1 SUBROUTINE DIIS_SIMPLE(LUEVEC,NVEC,NDIM,C) 2* 3* NVEC error-vectors are given on LUEVEC 4* Use DIIS to find optimal combination of these vectors 5* 6*. Simple and generic DIIS code - more simple than efficient 7*. All subspace matrices constructed without reuse 8* 9*. Jeppe Olsen, August 2005 10* 11 INCLUDE 'wrkspc.inc' 12*. Output : Optimial combination of directions 13 REAL*8 INPROD 14 DIMENSION C(NVEC) 15*. Non-redundant formalism : X = X(NVEC) + sum(i=1,NVEC-1)C(I)(X(I)-X(N)) 16* corresponding error E = E(NVEC) + sum(i=1,NVEC-1)C(I)(E(I)-E(N)) 17* 18* Minimizing norm of gradient corresponds to solving linear set of equations 19* 20* A C = RHS 21* A_ij = <E(I)-E(N)!E(J)-E(N)>, RHS_i = -<E(I)-E(N)!E(N)> 22* 23* 24 NTEST = 100 25* 26 IF(NTEST.GE.10) THEN 27 WRITE(6,*) ' Entering DIIS ' 28 WRITE(6,*) ' Number of vectors ', NVEC 29 END IF 30* 31 IDUM = 0 32 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'DIISS') 33*. A few scratch vectors 34 LEN = MAX(NDIM,NVEC-1) 35 CALL MEMMAN(KLVEC1,LEN,'ADDL ',2,'VC1_DI') 36 CALL MEMMAN(KLVEC2,NDIM,'ADDL ',2,'VC2_DI') 37 CALL MEMMAN(KLVEC3,NDIM,'ADDL ',2,'VC3_DI') 38*. Space for A and RHS and for inverting 39 CALL MEMMAN(KLA,(NVEC-1)**2,'ADDL ',2,'AMAT ') 40 CALL MEMMAN(KLRHS,(NVEC-1), 'ADDL ',2,'RHS ') 41 CALL MEMMAN(KLSCRM1,(NVEC-1)**2,'ADDL ',2,'SCRM1 ') 42 CALL MEMMAN(KLSCRM2,(NVEC-1)**2,'ADDL ',2,'SCRM2 ') 43 44*. Read E(NVEC) in VEC1 45 CALL REWINO(LUEVEC) 46 DO IVEC = 1, NVEC 47C? WRITE(6,*) ' Initial loop, IVEC = ', IVEC 48C VEC_FROM_DISC(VEC,LENGTH,IREW,LBLK,LU) 49 CALL VEC_FROM_DISC(WORK(KLVEC1),NDIM,0,-1,LUEVEC) 50 END DO 51* 52 ONE = 1.0D0 53 ONEM = -1.0D0 54 DO IVEC = 1, NVEC-1 55C? WRITE(6,*) ' IVEC = ', IVEC 56*. Read vector E(IVEC) in VEC2 57 CALL REWINO(LUEVEC) 58 DO KVEC = 1, IVEC 59 CALL VEC_FROM_DISC(WORK(KLVEC2),NDIM,0,-1,LUEVEC) 60 END DO 61*. E(IVEC)-E(NVEC) in VEC2 62 CALL VECSUM(WORK(KLVEC2),WORK(KLVEC2),WORK(KLVEC1), 63 & ONE,ONEM,NDIM) 64 WORK(KLRHS-1+IVEC) = INPROD(WORK(KLVEC2),WORK(KLVEC1),NDIM) 65 EEII = INPROD(WORK(KLVEC2),WORK(KLVEC2),NDIM) 66 WORK(KLA-1+(IVEC-1)*(NVEC-1)+IVEC) = EEII 67 DO JVEC = IVEC+1, NVEC-1 68C? WRITE(6,*) ' JVEC = ', JVEC 69*. Read vector E(JVEC) on VEC3 70 CALL VEC_FROM_DISC(WORK(KLVEC3),NDIM,0,-1,LUEVEC) 71*. E(JVEC) - E(NVEC) in VEC3 72 CALL VECSUM(WORK(KLVEC3),WORK(KLVEC3),WORK(KLVEC1), 73 & ONE,ONEM,NDIM) 74 EEIJ = INPROD(WORK(KLVEC2),WORK(KLVEC3),NDIM) 75 WORK(KLA-1+(IVEC-1)*(NVEC-1)+JVEC) = EEIJ 76 WORK(KLA-1+(JVEC-1)*(NVEC-1)+IVEC) = EEIJ 77 END DO 78 END DO 79*. Well, the RHS should contain a - so 80 CALL SCALVE(WORK(KLRHS),ONEM,NVEC-1) 81* 82 IF(NTEST.GE.100) THEN 83 WRITE(6,*) 'The <E(I)-E(N)!E(J)-E(N)> matrix in DIIS ' 84 CALL WRTMAT(WORK(KLA),NVEC-1,NVEC-1,NVEC-1,NVEC-1) 85 WRITE(6,*) 'The <E(I)-E(N)!E(N)> vector in DIIS ' 86 CALL WRTMAT(WORK(KLRHS),1,NVEC-1,1,NVEC-1) 87 END IF 88*. Solve linear equations 89*. Invert EE 90 IF(NVEC.GE.1) THEN 91 CALL INVERT_BY_DIAG2(WORK(KLA),WORK(KLSCRM1),WORK(KLSCRM2), 92 & WORK(KLVEC1),NVEC-1) 93 END IF 94 CALL MATVCB(WORK(KLA),WORK(KLRHS),C,NVEC-1,NVEC-1,0) 95*. IN C, we have coefficients for C(I),I.ne.NVEC, obtain C(NVEC) 96 C(NVEC) = 1.0D0-ELSUM(C,NVEC-1) 97*. Clean up : leave the file at end of vector NVEC 98 IF(NVEC.GE.2) THEN 99 CALL VEC_FROM_DISC(WORK(KLVEC2),NDIM,0,-1,LUEVEC) 100 END IF 101* 102 IF(NTEST.GE.10) THEN 103 WRITE(6,*) ' The C-coefficients obtained from DIIS ' 104 CALL WRTMAT(C,1,NVEC,1,NVEC) 105 END IF 106* 107C WRITE(6,*) ' Leaving DIIS ' 108 CALL MEMMAN(IDUM,IDUM,'FLUSM',IDUM,'DIISS') 109 RETURN 110 END 111 112 113 114c $Id$ 115