@@ -56,11 +56,11 @@ SUBROUTINE SETLAG(EOL)
5656!- ----------------------------------------------
5757! L o c a l V a r i a b l e s
5858!- ----------------------------------------------
59- INTEGER :: ITWICE, LIRAW, LI, LIP1, NAKLI , LJRAW, LJ, IECCLI, L1, L2, &
59+ INTEGER :: ITWICE, LIRAW, LI, LIP1, NAKLJ , LJRAW, LJ, IECCLI, L1, L2, &
6060 JLAST, MLAST, M, J, I
6161 REAL (DOUBLE), DIMENSION (NNNP) :: YPJ, YPM, XPJ, XPM, XQJ, XQM
6262 REAL (DOUBLE) :: EPS, UCFJ, UCFM, RESULT, RIJM, QDIF, OBQDIF, OBQSUM
63- LOGICAL :: FIRST, FIXLI, FIXLJ, FULLI, FULLJ
63+ LOGICAL :: FIRST, FIXLI, FIXLJ, FULLI, FULLJ, VLI, VLJ
6464!- ----------------------------------------------
6565!
6666 DATA FIRST/ .TRUE. /
@@ -85,65 +85,28 @@ SUBROUTINE SETLAG(EOL)
8585 EPS = ACCY* 0.01D0 ! criterion to see if an orb is occupied
8686 DO ITWICE = 1 , 2
8787 NEC = 0
88- ! IF (ITWICE /= 2) THEN
89- ! DO LIRAW = 1, NW - 1
90- ! LI = IORDER(LIRAW)
91- ! LIP1 = MAX(NCORE,LIRAW) + 1
92- ! NAKLI = NAK(LI)
93- ! FIXLI = LFIX(LI)
94- ! FULLI = ABS(UCF(LI)-DBLE(NKJ(LI)+1)) < EPS
95- ! DO LJRAW = LIP1, NW
96- ! LJ = IORDER(LJRAW)
97- ! FIXLJ = LFIX(LJ)
98- ! FULLJ = ABS(UCF(LJ)-DBLE(NKJ(LJ)+1)) < EPS
99- ! IF (.NOT.(NAK(LJ)==NAKLI .AND. .NOT.(FIXLI .AND. FIXLJ)&
100- ! .AND. .NOT.(FULLI .AND. FULLJ))) CYCLE
101- ! NEC = NEC + 1
102- ! CYCLE
103- ! !*** Encode index at 2nd round ***
104- ! END DO
105- ! END DO
106- ! ELSE
107- ! DO LIRAW = 1, NW - 1
108- ! LI = IORDER(LIRAW)
109- ! LIP1 = MAX(NCORE,LIRAW) + 1
110- ! NAKLI = NAK(LI)
111- ! FIXLI = LFIX(LI)
112- ! FULLI = ABS(UCF(LI)-DBLE(NKJ(LI)+1)) < EPS
113- ! !*** Encode index at 2nd round ***
114- ! DO LJRAW = LIP1, NW
115- ! LJ = IORDER(LJRAW)
116- ! FIXLJ = LFIX(LJ)
117- ! FULLJ = ABS(UCF(LJ)-DBLE(NKJ(LJ)+1)) < EPS
118- ! IF (.NOT.(NAK(LJ)==NAKLI .AND. .NOT.(FIXLI .AND. FIXLJ)&
119- ! .AND. .NOT.(FULLI .AND. FULLJ))) CYCLE
120- ! NEC = NEC + 1
121- ! !*** Encode index at 2nd round ***
122- ! IECC(NEC) = LI + KEY*LJ
123- ! END DO
124- ! END DO
125- ! ENDIF
126- DO LIraw = 1 , NW - 1
127- LI = iorder(LIraw)
128- LIP1 = MAX (NCORE, LIraw) + 1
129- NAKLI = NAK(LI)
130- FIXLI = LFIX(LI)
131- FULLI = ABS ( UCF(LI)- DBLE (NKJ(LI)+ 1 ) ) .LT. EPS
132- DO LJraw = LIP1, NW
88+ DO LJraw = NCORE+1 , NW
13389 LJ = iorder(LJraw)
134- FIXLJ = LFIX(LJ)
90+ NAKLJ = NAK(LJ)
91+ VLJ = .NOT. LFIX(LJ)
13592 FULLJ = ABS ( UCF(LJ)- DBLE (NKJ(LJ)+ 1 ) ) .LT. EPS
136- IF ( (NAK(LJ) .EQ. NAKLI) .AND. &
137- (.NOT. (FIXLI .AND. FIXLJ)) .AND. &
138- (.NOT. (FULLI .AND. FULLJ)) ) THEN
93+ DO LIraw = 1 , LJraw-1
94+ LI = iorder(LIraw)
95+ VLI = .NOT. LFIX(LI)
96+ FULLI = ABS ( UCF(LI)- DBLE (NKJ(LI)+ 1 ) ) .LT. EPS
97+ IF (NAK(LI) .EQ. NAKLJ) then
98+ If (VLI .OR. VLJ ) then ! at least one orbital varied
99+ ! but not both (varied and full)
100+ If (.NOT. ((VLI .AND. VLJ) .AND. (FULLI .AND. FULLJ))) THEN
139101 NEC = NEC + 1
140102 ! *** Encode index at 2nd round ***
141103 IF (itwice == 2 ) IECC(NEC) = LI + KEY * LJ
142104 ENDIF
105+ ENDIF
106+ ENDIF
143107 ENDDO
144108 ENDDO
145109
146-
147110 IF (ITWICE== 1 .AND. NEC> 0 ) THEN
148111 CALL ALLOC (ECV, NEC, ' ECV' , ' SETLAG' )
149112 CALL ALLOC (IECC, NEC, ' IECC' , ' SETLAG' )
@@ -173,7 +136,6 @@ SUBROUTINE SETLAG(EOL)
173136 FIRST = .FALSE.
174137 ENDIF
175138
176- ! FF+GG 12/07/05
177139! Lagrange multipliers need to be computed also on the first call
178140! RETURN
179141
@@ -217,11 +179,12 @@ SUBROUTINE SETLAG(EOL)
217179 IF (LFIX(M)) THEN
218180 TA(1 ) = 0.D0
219181 DO I = 2 , MTP
220- TA(I) = RPOR(I)* ((PF(I,M)* XQJ(I)- QF(I,M)* XPJ(I))* C+ (PF(I,M) * PF(I &
221- , J)+ QF(I,M)* QF(I,J))* YPJ(I))
182+ TA(I) = RPOR(I)* ((PF(I,M)* XQJ(I)- QF(I,M)* XPJ(I))* C+ &
183+ (PF(I,M) * PF(I, J)+ QF(I,M)* QF(I,J))* YPJ(I))
222184 END DO
223185
224186 CALL QUAD (RESULT)
187+ ! note ..rinti(m,j,1) is symmetric in (j,m)
225188 RIJM = RINTI(M,J,1 )
226189 ECV(LI) = (RESULT - RIJM)* UCFJ
227190
@@ -236,8 +199,8 @@ SUBROUTINE SETLAG(EOL)
236199 ELSE IF (LFIX(J)) THEN
237200 TA(1 ) = 0.D0
238201 DO I = 2 , MTP
239- TA(I) = RPOR(I)* ((PF(I,J)* XQM(I)- QF(I,J)* XPM(I))* C+ (PF(I,J) * PF(I &
240- , M)+ QF(I,J)* QF(I,M))* YPM(I))
202+ TA(I) = RPOR(I)* ((PF(I,J)* XQM(I)- QF(I,J)* XPM(I))* C+ &
203+ (PF(I,J) * PF(I, M)+ QF(I,J)* QF(I,M))* YPM(I))
241204 END DO
242205
243206! start dbg
@@ -258,20 +221,22 @@ SUBROUTINE SETLAG(EOL)
258221! WRITE (81,*)RESULT, RIJM, UCFJ, ECV, r(i), rp(i)
259222! end dbg
260223
261-
262224 ELSE
225+
263226 QDIF = ABS ((UCFJ - UCFM)/ MAX (UCFJ,UCFM))
227+
264228 IF (QDIF > P001) THEN
265229 OBQDIF = 1.D0 / UCFJ - 1.D0 / UCFM
266230 TA(1 ) = 0.D0
267231 DO I = 2 , MTP
268- TA(I) = RPOR(I)* ((PF(I,M)* XQJ(I)- QF(I,M)* XPJ(I)- PF(I,J) * XQM(I &
269- ) + QF (I,J)* XPM (I)) * C + (YPJ(I) - YPM(I)) * (PF(I,M) * PF(I,J) + QF(I, &
270- M)* QF(I,J)))
232+ TA(I) = RPOR(I)* ((PF(I,M)* XQJ(I)- QF(I,M)* XPJ(I) &
233+ - PF (I,J)* XQM (I)+ QF(I,J) * XPM(I)) * C &
234+ + (YPJ(I) - YPM(I)) * (PF(I, M)* PF(I,J) + QF(I,M) * QF(I, J)))
271235 END DO
272236
273237 CALL QUAD (RESULT)
274238 ECV(LI) = RESULT/ OBQDIF
239+
275240! start dbg
276241! WRITE (81,*)'3, RESULT, OBQDIF, ECV, TA'
277242! WRITE (81,*)RESULT, OBQDIF, ECV
@@ -280,19 +245,22 @@ SUBROUTINE SETLAG(EOL)
280245! ENDDO
281246! end dbg
282247
283-
284248 ELSE
249+
285250 OBQSUM = 1.D0 / UCFJ + 1.D0 / UCFM
251+
286252 TA(1 ) = 0.D0
287253 DO I = 2 , MTP
288- TA(I) = RPOR(I)* ((PF(I,M)* XQJ(I)- QF(I,M)* XPJ(I)+ PF(I,J) * XQM(I &
289- ) - QF (I,J)* XPM (I)) * C + (YPJ(I) + YPM(I)) * (PF(I,M) * PF(I,J) + QF(I, &
290- M)* QF(I,J)))
254+ TA(I) = RPOR(I)* ((PF(I,M)* XQJ(I)- QF(I,M)* XPJ(I) &
255+ + PF (I,J)* XQM (I)- QF(I,J) * XPM(I)) * C &
256+ + (YPJ(I) + YPM(I)) * (PF(I, M)* PF(I,J) + QF(I,M) * QF(I, J)))
291257 END DO
292258
293259 CALL QUAD (RESULT)
294- RIJM = RINTI(M,J,1 ) ! / nprocs
295- ECV(LI) = (RESULT - 2.D0 * RIJM)/ OBQSUM
260+ RIJM = RINTI(M,J,1 )
261+ ECV(LI) = (RESULT - 2.0 * RIJM)/ OBQSUM
262+
263+
296264! start dbg
297265! WRITE (81,*)'4, RESULT, RIUJM, OBQSUM, ECV, TA'
298266! WRITE (81,*)RESULT, RIUJM, OBQSUM, ECV
@@ -308,10 +276,6 @@ SUBROUTINE SETLAG(EOL)
308276
309277 END DO
310278
311- ! db close(81)
312- ! db close(82)
313-
314-
315279 302 FORMAT (/ ,' Lagrange multipliers are not required' )
316280 304 FORMAT (/ ,' Include Lagrange multipliers between:' / )
317281 305 FORMAT (13X ,2 (2X ,1I2 ,1A2 ))
0 commit comments