hdiff output

r29119/finalio.f90 2015-11-17 23:33:06.548592573 +0000 r29118/finalio.f90 2015-11-17 23:33:07.916610920 +0000
587:                  DO J5=1,3587:                  DO J5=1,3
588:                     P(J5) = QMINP(J1,3*(J4-1)+J5)                  588:                     P(J5) = QMINP(J1,3*(J4-1)+J5)                  
589:                     IF(PERIODIC) THEN ! wrap back into the box589:                     IF(PERIODIC) THEN ! wrap back into the box
590:                        P(J5)=P(J5)-BOX3D(J5)*ANINT(P(J5)/BOX3D(J5))590:                        P(J5)=P(J5)-BOX3D(J5)*ANINT(P(J5)/BOX3D(J5))
591:                     ENDIF591:                     ENDIF
592:                  ENDDO592:                  ENDDO
593:                  WRITE(MYUNIT2,'(A,3(1X,F20.10))',ADVANCE='NO') &593:                  WRITE(MYUNIT2,'(A,3(1X,F20.10))',ADVANCE='NO') &
594:                       ATOM_TYPE,(P(J5), J5=1,3)594:                       ATOM_TYPE,(P(J5), J5=1,3)
595:                  IF(STRESST) THEN ! Tack on local stresses595:                  IF(STRESST) THEN ! Tack on local stresses
596:                     IF(STRESS_MODE==2) THEN596:                     IF(STRESS_MODE==2) THEN
597:                        ! Just print the uniqe elements of the stress597:                        ! Just pring the uniqe elements of the stress
598:                        ! tensor598:                        ! tensor
599:                        WRITE(MYUNIT2,'(6(1X,F10.5))',ADVANCE='YES') &599:                        WRITE(MYUNIT2,'(6(1X,F10.5))',ADVANCE='YES') &
600:                             STRESS(J4,1,1),STRESS(J4,2,2),STRESS(J4,3,3),&600:                             STRESS(J4,1,1),STRESS(J4,2,2),STRESS(J4,3,3),&
601:                             STRESS(J4,1,2),STRESS(J4,1,3),STRESS(J4,2,3)601:                             STRESS(J4,1,2),STRESS(J4,1,3),STRESS(J4,2,3)
602:                     ELSEIF(STRESS_MODE==1) THEN602:                     ELSEIF(STRESS_MODE==1) THEN
603:                        WRITE(MYUNIT2,'(2(1X,F10.5))',ADVANCE='YES') &603:                        WRITE(MYUNIT2,'(2(1X,F10.5))',ADVANCE='YES') &
604:                             -(STRESS(J4,1,1)+STRESS(J4,2,2)+STRESS(J4,3,3))/3.0D0, &604:                             -(STRESS(J4,1,1)+STRESS(J4,2,2)+STRESS(J4,3,3))/3.0D0, &
605:                             DSQRT( (STRESS(J4,1,1)-STRESS(J4,2,2))**2 + &605:                             DSQRT( (STRESS(J4,1,1)-STRESS(J4,2,2))**2 + &
606:                             (STRESS(J4,2,2)-STRESS(J4,3,3))**2 + &606:                             (STRESS(J4,2,2)-STRESS(J4,3,3))**2 + &
607:                             (STRESS(J4,3,3)-STRESS(J4,1,1))**2 )/3.0D0607:                             (STRESS(J4,3,3)-STRESS(J4,1,1))**2 )/3.0D0


r29119/MGupta.f90 2015-11-17 23:33:06.360590050 +0000 r29118/MGupta.f90 2015-11-17 23:33:07.728608396 +0000
 15: ! 15: !
 16: !   You should have received a copy of the GNU General Public License 16: !   You should have received a copy of the GNU General Public License
 17: !   along with this program; IF not, write to the Free Software 17: !   along with this program; IF not, write to the Free Software
 18: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 18: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 19: ! 19: !
 20: !****************************************************************** 20: !******************************************************************
 21: !   Multi-component Gupta potential. 21: !   Multi-component Gupta potential.
 22: !   -ds656 (Dec 2014) 22: !   -ds656 (Dec 2014)
 23: !****************************************************************** 23: !******************************************************************
 24: ! 24: !
 25: SUBROUTINE MGUPTA(X,GRAD,POT,GRADT,HESST,STRESST) 25: SUBROUTINE MGUPTA (X,V,PG,GRADT,STRESST)
 26:   ! 26:   !
 27:   USE COMMONS, ONLY : NATOMS, VT, NSPECIES, ATOMLISTS,& 27:   USE COMMONS, ONLY : NATOMS,VT,ATOMLISTS,NSPECIES,STRESS
 28:        INVATOMLISTS, STRESS 
 29:   USE MODHESS, ONLY : HESS 
 30:   USE POT_PARAMS, ONLY : MGUPTA_A, MGUPTA_XI, MGUPTA_P, & 28:   USE POT_PARAMS, ONLY : MGUPTA_A, MGUPTA_XI, MGUPTA_P, &
 31:        MGUPTA_Q, MGUPTA_R0, MGUPTA_M2Q, MGUPTA_XI_SQ, & 29:        MGUPTA_Q, MGUPTA_R0
 32:        MGUPTA_MP_DIVBY_R0, MGUPTA_M2Q_DIVBY_R0 
 33:   ! 30:   !
 34:   IMPLICIT NONE 31:   IMPLICIT NONE
 35:   ! 32:   !
 36:   ! Passed variables ----------------------------------------- 33:   ! Passed variables -----------------------------------------
 37:   LOGICAL, INTENT(IN) :: GRADT, HESST, STRESST 34:   LOGICAL, INTENT(IN) :: GRADT, STRESST
 38:   DOUBLE PRECISION, INTENT(IN) :: X(3*NATOMS) 35:   DOUBLE PRECISION, INTENT(IN) :: X(3*NATOMS)
 39:   DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), POT 36:   DOUBLE PRECISION, INTENT(OUT) :: V(3*NATOMS), PG
 40:   ! 37:   !
 41:   ! Local variables ----------------------------------------- 38:   ! Local variables -----------------------------------------
 42:   LOGICAL :: SELF 39:   LOGICAL :: SELF
 43:   INTEGER :: I,J,J1,J2,J3,J13,J23,J33,T1,T2,T3,G1,G2,G3,& 40:   INTEGER :: I,J,J1,J2,J13,J23,ATOMTYPE,MBI,T1,T2,G1,G2,K1,K2,&
 44:        K1,K2,K3,LI1,UI1,LI2,UI2,G2START 41:        LI1,UI1,LI2,UI2,G2START
 45:   DOUBLE PRECISION :: DIST, DX(3), RHO(NATOMS), DUMMY, & 42:   DOUBLE PRECISION :: DIST, DX(3), &
 46:        RIJ(NATOMS,NATOMS), DUDR, D2UDR2, IGRAD(3), & 43:        GRHO(NATOMS), DUDR, DUMMY, DISTANCE_MATRIX(NATOMS,NATOMS), &
 47:        REP_IJ(NATOMS,NATOMS), ATT_IJ(NATOMS,NATOMS), & 44:        VTEMP(NATOMS,3), PWR1, PWR2, GRAD(3)
 48:        REP_DUM1, REP_DUM2, REP_DUM3, ATT_DUM1, ATT_DUM2, & 
 49:        ATT_DUM3, RHO_DUM, RHO4HESS(NATOMS) 
 50:   ! 45:   !
 51:   ! Do some initialisations 46:   ! Do some initialisations
 52:   RIJ(1:NATOMS,1:NATOMS)=0.0D0 47:   DISTANCE_MATRIX(1:NATOMS,1:NATOMS)=0.0D0
 53:   REP_IJ(1:NATOMS,1:NATOMS)=0.0D0 48:   GRHO(1:NATOMS)=0.0D0
 54:   ATT_IJ(1:NATOMS,1:NATOMS)=0.0D0 
 55:   RHO(1:NATOMS)=0.0D0 
 56:   RHO4HESS(1:NATOMS)=0.0D0 
 57:   VT(1:NATOMS) = 0.0D0 49:   VT(1:NATOMS) = 0.0D0
 58:   ! 50:   !
  51:   MBI = 2 ! Many-Body Index (2 for Gupta; 1 for Morse-like pair pot.)
  52:   PWR1 = 1.0D0/DBLE(MBI)
  53:   PWR2 = PWR1-1.0D0
  54:   !
 59:   ! Double loop over atom types 55:   ! Double loop over atom types
 60:   LT1: DO T1=1,NSPECIES(0) 56:   LT1: DO T1=1,NSPECIES(0)
 61:      LT2: DO T2=1,T1 ! T1,NSPECIES(0) 57:      LT2: DO T2=1,T1 ! T1,NSPECIES(0)
 62:         ! 58:         !
 63:         ! Double loop over mobile and frozen atom groups 59:         ! Double loop over mobile and frozen atom groups
 64:         LG1: DO G1 = 1,2 ! 1 -> mobile; 2 -> frozen 60:         LG1: DO G1 = 1,2 ! 1 -> mobile; 2 -> frozen
 65:            ! 61:            !
 66:            IF(T1 == T2) THEN 62:            IF(T1 == T2) THEN
 67:               G2START = G1 63:               G2START = G1
 68:            ELSE 64:            ELSE
103:                     ! 99:                     !
104:                     !TESTCOUNT=TESTCOUNT+1100:                     !TESTCOUNT=TESTCOUNT+1
105:                     !WRITE(*,*) "BGupta> TEST:", J1,J2,T1,T2101:                     !WRITE(*,*) "BGupta> TEST:", J1,J2,T1,T2
106:                     !102:                     !
107:                     DIST=0.0D0103:                     DIST=0.0D0
108:                     DO I=1,3104:                     DO I=1,3
109:                        DIST = DIST + (X(J23+I)-X(J13+I))**2105:                        DIST = DIST + (X(J23+I)-X(J13+I))**2
110:                     ENDDO106:                     ENDDO
111:                     DIST = DSQRT(DIST)107:                     DIST = DSQRT(DIST)
112:                     !108:                     !
113:                     RIJ(J2,J1)=DIST ! store distance109:                     DISTANCE_MATRIX(J2,J1)=DIST ! store distance
114:                     RIJ(J1,J2)=DIST ! impose symmetry110:                     DISTANCE_MATRIX(J1,J2)=DIST ! impose symmetry
115:                     !111:                     !
116:                     DUMMY = MGUPTA_A(T1,T2)*&112:                     DIST=DIST/MGUPTA_R0(T1,T2)
117:                          DEXP(MGUPTA_MP_DIVBY_R0(T1,T2)*DIST + & 
118:                          MGUPTA_P(T1,T2)) 
119:                     REP_IJ(J1,J2) = DUMMY 
120:                     REP_IJ(J2,J1) = DUMMY 
121:                     !113:                     !
 114:                     DUMMY = MGUPTA_A(T1,T2)*DEXP(MGUPTA_P(T1,T2)*(1-DIST))
122:                     VT(J1) = VT(J1) + DUMMY115:                     VT(J1) = VT(J1) + DUMMY
123:                     VT(J2) = VT(J2) + DUMMY116:                     VT(J2) = VT(J2) + DUMMY
124:                     !117:                     !
125:                     DUMMY = MGUPTA_XI_SQ(T1,T2)*&118:                     ! calculate many-body potential term:
126:                          DEXP(MGUPTA_M2Q_DIVBY_R0(T1,T2)*DIST - &119:                     DUMMY = DEXP(2.0D0*MGUPTA_Q(T1,T2)*(1-DIST)) 
127:                          MGUPTA_M2Q(T1,T2))120:                     DUMMY = DUMMY * MGUPTA_XI(T1,T2)**MBI
128:                     ATT_IJ(J1,J2) = DUMMY121:                     GRHO(J1) = GRHO(J1) + DUMMY
129:                     ATT_IJ(J2,J1) = DUMMY122:                     GRHO(J2) = GRHO(J2) + DUMMY
130:                     ! 
131:                     ! Accumulate many-body potential term: 
132:                     RHO(J1) = RHO(J1) + DUMMY 
133:                     RHO(J2) = RHO(J2) + DUMMY 
134:                     ! 
135:                     IF (GRADT.OR.HESST.OR.STRESST) THEN 
136:                        ! 
137:                        ! Precompute derivatives of rep and att here 
138:                        REP_IJ(J1,J2) = MGUPTA_MP_DIVBY_R0(T1,T2)*& 
139:                             REP_IJ(J1,J2)*2.0D0  ! need the 2 here! 
140:                        REP_IJ(J2,J1) = REP_IJ(J1,J2) 
141:                        ATT_IJ(J1,J2) = MGUPTA_M2Q_DIVBY_R0(T1,T2)*& 
142:                             ATT_IJ(J1,J2) 
143:                        ATT_IJ(J2,J1) = ATT_IJ(J1,J2) 
144:                        ! 
145:                     ENDIF 
146:                     !123:                     !
147:                  ENDDO LA2 ! Double loop over atoms124:                  ENDDO LA2 ! Double loop over atoms
148:               ENDDO LA1125:               ENDDO LA1
149:               !126:               !
150:            ENDDO LG2 ! Span dynamic (1) and frozen (2) groups only127:            ENDDO LG2 ! Span dynamic (1) and frozen (2) groups only
151:         ENDDO LG1128:         ENDDO LG1
152:         !129:         !
153:      ENDDO LT2 ! Double loop over atom types130:      ENDDO LT2 ! Double loop over atom types
154:   ENDDO LT1131:   ENDDO LT1
155:   !132:   !
156:   ! Now store per-atom energies in array VT and sum over all atoms133:   ! Now store per-atom energies in array VT and sum over all atoms
157:   !134:   !
158:   POT=0.0D0135:   PG=0.0D0
159:   DO G1=1,2 ! Groups 1 and 2136:   DO G1=1,2 ! Groups 1 and 2
160:      DO T1=1,NSPECIES(0) ! Atom types137:      DO T1=1,NSPECIES(0) ! Atom types
161:         DO J2=1,ATOMLISTS(T1,G1,0) 138:         DO J2=1,ATOMLISTS(T1,G1,0) 
162:            J1=ATOMLISTS(T1,G1,J2)139:            J1=ATOMLISTS(T1,G1,J2)
163:            RHO(J1)=DSQRT(RHO(J1)) ! square root density140:            !GRHO(J1)=DSQRT(GRHO(J1)) ! square root density
164:            VT(J1) = VT(J1) - RHO(J1) ! accummulate per-atom energy141:            VT(J1) = VT(J1) - GRHO(J1)**PWR1 ! accummulate per-atom energy
165:            POT = POT + VT(J1)       ! accumulate potential energy142:            PG = PG + VT(J1)       ! accumulate potential energy
166:            IF(GRADT.OR.GRADT.OR.STRESST) THEN 
167:               RHO(J1)=1.0D0/RHO(J1)  
168:               IF(HESST) RHO4HESS(J1) = 0.25D0*RHO(J1)**3 
169:               RHO(J1)=0.5D0*RHO(J1) 
170:            ENDIF 
171:         ENDDO143:         ENDDO
172:      ENDDO144:      ENDDO
173:   ENDDO145:   ENDDO
174:   !146:   !
175:   ! Calculate gradient terms, if required147:   ! Calculate gradient terms, if required
176:   !148:   !
177:   IF (GRADT.OR.HESST.OR.STRESST) THEN149:   IF (GRADT) THEN
178:      ! 150:      ! 
179:      GRAD(:) = 0.0D0151:      V(:) = 0.0d0
180:      IF(HESST) HESS(:,:) = 0.0D0152:      VTEMP(:,:)=0.0D0
181:      IF(STRESST) STRESS(:,:,:) = 0.0D0 
182:      !153:      !
183:      ! Double loop over atom types154:      ! Double loop over atom types
184:      LT3: DO T1=1,NSPECIES(0)155:      LT3: DO T1=1,NSPECIES(0)
185:         LT4: DO T2=1,T1 156:         LT4: DO T2=1,T1 
186:            !157:            !
187:            ! Double loop over free / frozen groups158:            ! Double loop over free / frozen groups
188:            LG3: DO G1 = 1,2 ! loop over goups 1 and 2159:            LG3: DO G1 = 1,2 ! loop over goups 1 and 2
189:               !160:               !
190:               IF(T1 == T2) THEN161:               IF(T1 == T2) THEN
191:                  G2START = G1162:                  G2START = G1
203:                     UI1 = UI1 - 1 ! Exclude self-interaction174:                     UI1 = UI1 - 1 ! Exclude self-interaction
204:                  ELSE175:                  ELSE
205:                     SELF = .FALSE.176:                     SELF = .FALSE.
206:                  ENDIF177:                  ENDIF
207:                  !178:                  !
208:                  LA3: DO K1=LI1,UI1179:                  LA3: DO K1=LI1,UI1
209:                     !180:                     !
210:                     J1 = ATOMLISTS(T1,G1,K1) ! actual atom 1 index181:                     J1 = ATOMLISTS(T1,G1,K1) ! actual atom 1 index
211:                     J13 = 3*(J1-1)182:                     J13 = 3*(J1-1)
212:                     !183:                     !
 184:                     ! store reciprocal of RHO for atom J1
 185:                     DUMMY=GRHO(J1)**PWR2 
 186:                     !
213:                     ! 2nd atom loop bounds depend on group.187:                     ! 2nd atom loop bounds depend on group.
214:                     ! Can this IF block be moved outside of loop A1?188:                     ! Can this IF block be moved outside of loop A1?
215:                     IF(SELF) THEN189:                     IF(SELF) THEN
216:                        LI2=K1+1 190:                        LI2=K1+1 
217:                        UI2=UI1+1191:                        UI2=UI1+1
218:                     ELSE192:                     ELSE
219:                        LI2=1193:                        LI2=1
220:                        UI2=ATOMLISTS(T2,G2,0)194:                        UI2=ATOMLISTS(T2,G2,0)
221:                     ENDIF195:                     ENDIF
222:                     !196:                     !
223:                     LA4: DO K2=LI2,UI2197:                     LA4: DO K2=LI2,UI2
224:                        !198:                        !
225:                        J2 = ATOMLISTS(T2,G2,K2) ! actual atom 2 index199:                        J2 = ATOMLISTS(T2,G2,K2) ! actual atom 2 index
226:                        J23 = 3*(J2-1)200:                        J23 = 3*(J2-1)
227:                        !201:                        !
228:                        RHO_DUM = RHO(J1) + RHO(J2) ! factor of 1/2 is included202:                        DIST=DISTANCE_MATRIX(J1,J2) ! recall distance 
229:                        DUDR = REP_IJ(J1,J2) - RHO_DUM*ATT_IJ(J1,J2)203:                        DIST=DIST/MGUPTA_R0(T1,T2)
230:                        ! Want 1/R*dU/dR204:                        !
231:                        DIST = RIJ(J1,J2)205:                        ! DUDR is basically 1/R*dU/dR
232:                        DUDR = DUDR/DIST206:                        DUDR = 2.0D0*( PWR1*MGUPTA_Q(T1,T2)* &
 207:                             (MGUPTA_XI(T1,T2)**MBI)* &
 208:                             (DUMMY+GRHO(J2)**PWR2)* &
 209:                             DEXP(2.0D0*MGUPTA_Q(T1,T2)*(1-DIST)) - &
 210:                             MGUPTA_A(T1,T2)*MGUPTA_P(T1,T2)* &
 211:                             DEXP(MGUPTA_P(T1,T2)*(1-DIST))) &
 212:                             /(MGUPTA_R0(T1,T2)**2*DIST) 
233:                        !213:                        !
234:                        ! Cartesian components of distance214:                        ! Cartesian components of distance
235:                        DO I=1,3215:                        DO I=1,3
236:                           DX(I)=X(J13+I)-X(J23+I)216:                           DX(I)=X(J13+I)-X(J23+I)
237:                           IGRAD(I)=DUDR*DX(I) 217:                           GRAD(I)=DUDR*DX(I) 
238:                           GRAD(J13+I) = GRAD(J13+I) + IGRAD(I)218:                           VTEMP(J1,I)=VTEMP(J1,I)+GRAD(I)
239:                           GRAD(J23+I) = GRAD(J23+I) - IGRAD(I)219:                           ! impose symmetry on gradient
 220:                           VTEMP(J2,I)=VTEMP(J2,I)-GRAD(I) 
240:                        ENDDO221:                        ENDDO
241:                        !222:                        !
242:                        IF(STRESST) THEN ! Accumulate local stresses223:                        IF(STRESST) THEN ! Accumulate local stresses
243:                           DO I=1,3224:                           DO I=1,3
244:                              DO J=I,3225:                              DO J=I,3
245:                                 DUMMY = IGRAD(I)*DX(J) 
246:                                 STRESS(J1,I,J) = STRESS(J1,I,J) + &226:                                 STRESS(J1,I,J) = STRESS(J1,I,J) + &
247:                                      DUMMY227:                                      GRAD(I)*DX(J)
 228:                                 ! + or - ?!? I think +
248:                                 STRESS(J2,I,J) = STRESS(J2,I,J) + &229:                                 STRESS(J2,I,J) = STRESS(J2,I,J) + &
249:                                      DUMMY230:                                      GRAD(I)*DX(J)
250:                                 STRESS(0,I,J) = STRESS(0,I,J) + & 
251:                                      DUMMY 
252:                              ENDDO231:                              ENDDO
253:                           ENDDO232:                           ENDDO
254:                        ENDIF233:                        ENDIF
255:                        !234:                        !
256:                        IF(HESST) THEN 
257:                           ! 
258:                           DUMMY=1.0D0/DIST**2 ! Reset DUMMY 
259:                           REP_DUM1=REP_IJ(J1,J2)/DIST 
260:                           REP_DUM2=MGUPTA_MP_DIVBY_R0(T1,T2)*& 
261:                                REP_IJ(J1,J2) 
262:                           REP_DUM2=(REP_DUM2-REP_DUM1)*DUMMY 
263:                           ATT_DUM1=ATT_IJ(J1,J2)/DIST 
264:                           ATT_DUM2=MGUPTA_M2Q_DIVBY_R0(T1,T2)*& 
265:                                ATT_IJ(J1,J2) 
266:                           ATT_DUM2=(ATT_DUM2-ATT_DUM1)*DUMMY 
267:                           ! 
268:                           DO I = 1,3 ! coordinates of atom J1 
269:                              DO J= 1,3 ! coordinate of atom J2 (/=J1 by construction!) 
270:                                 !                                 
271:                                 DUMMY = DX(I)*DX(J) ! Reset DUMMY 
272:                                 ! 
273:                                 REP_DUM3 = REP_DUM2*DUMMY 
274:                                 ATT_DUM3 = ATT_DUM2*DUMMY 
275:                                 ! 
276:                                 IF (I == J) THEN 
277:                                    REP_DUM3 = REP_DUM3 + REP_DUM1 
278:                                    ATT_DUM3 = ATT_DUM3 + ATT_DUM1 
279:                                 ENDIF 
280:                                 ! 
281:                                 D2UDR2 = REP_DUM3 - RHO_DUM*ATT_DUM3 
282:                                 DUDR = D2UDR2 
283:                                 ! 
284:                                 ! Now go through triplets... FECK! 
285:                                 DO T3=1,NSPECIES(0) 
286:                                    DO G3 = 1,2 
287:                                       DO K3=1,ATOMLISTS(T3,G3,0)  
288:                                          ! 
289:                                          J3 = ATOMLISTS(T3,G3,K3) 
290:                                          J33 = 3*(J3-1) 
291:                                          ! 
292:                                          IF ((J3.NE.J2).AND.(J3.NE.J1)) THEN                                             
293:                                             DUMMY = RHO4HESS(J3)*ATT_IJ(J1,J3)*ATT_IJ(J2,J3)/& 
294:                                                  (RIJ(J1,J3)*RIJ(J2,J3))                                             
295:                                             D2UDR2 = D2UDR2 - & 
296:                                                  DUMMY*(X(J13+I)-X(J33+I))*(X(J23+J)-X(J33+J))                                             
297:                                             DUDR = DUDR - & 
298:                                                  DUMMY*(X(J23+I)-X(J33+I))*(X(J13+J)-X(J33+J)) 
299:                                          ENDIF 
300:                                          ! 
301:                                          IF (J3.NE.J2) THEN 
302:                                             DUMMY = RHO4HESS(J2)* & 
303:                                                  ATT_IJ(J1,J2)*ATT_IJ(J3,J2)/& 
304:                                                  (RIJ(J1,J2)*RIJ(J3,J2)) 
305:                                             D2UDR2 = D2UDR2 + & 
306:                                                  DUMMY*(X(J13+I)-X(J23+I))*(X(J33+J)-X(J23+J)) 
307:                                             DUDR = DUDR + & 
308:                                                  DUMMY*(X(J33+I)-X(J23+I))*(X(J13+J)-X(J23+J)) 
309:                                          ENDIF 
310:                                          ! 
311:                                          IF (J3.NE.J1) THEN 
312:                                             DUMMY = RHO4HESS(J1)* & 
313:                                                  ATT_IJ(J3,J1)*ATT_IJ(J2,J1)/&                                                  
314:                                                  (RIJ(J3,J1)*RIJ(J2,J1)) 
315:                                             D2UDR2 = D2UDR2 + & 
316:                                                  DUMMY*(X(J33+I)-X(J13+I))*(X(J23+J)-X(J13+J)) 
317:                                             DUDR = DUDR + & 
318:                                                  DUMMY*(X(J23+I)-X(J13+I))*(X(J33+J)-X(J13+J)) 
319:                                          ENDIF 
320:                                          ! 
321:                                       ENDDO 
322:                                    ENDDO 
323:                                 ENDDO 
324:                                 ! 
325:                                 ! Accumulate diagonal blocks first 
326:                                 HESS(J13+I, J13+J) = & 
327:                                      HESS(J13+I, J13+J) + D2UDR2 
328:                                 HESS(J23+I, J23+J) = & 
329:                                      HESS(J23+I, J23+J) + DUDR ! Note D2UDR2 /= DUDR 
330:                                 ! 
331:                                 ! Now compute (not accumulate!) off-diagonal blocks, 
332:                                 ! Imposing symmetry over J1 and J2. 
333:                                 D2UDR2 = -D2UDR2 
334:                                 HESS(J13+I, J23+J) = D2UDR2 
335:                                 HESS(J23+J, J13+I) = D2UDR2 
336:                                 ! 
337:                              ENDDO 
338:                           ENDDO 
339:                           ! 
340:                        ENDIF ! End of Hessian(J1,J2) block 
341:                        ! 
342:                     ENDDO LA4235:                     ENDDO LA4
343:                  ENDDO LA3236:                  ENDDO LA3
344:                  !237:                  !
345:               ENDDO LG4238:               ENDDO LG4
346:            ENDDO LG3239:            ENDDO LG3
347:            !240:            !
348:         ENDDO LT4241:         ENDDO LT4
349:      ENDDO LT3242:      ENDDO LT3
350:      !243:      !
351:      ! Finally, zero gradiend for group 2 atoms244:      ! Finally, sum the gradient for group 1 atoms only
352:      DO T1=1,NSPECIES(0) ! Atom types245:      DO T1=1,NSPECIES(0) ! Atom types
353:         DO J2=1,ATOMLISTS(T1,2,0) ! Group 2246:         DO J2=1,ATOMLISTS(T1,1,0) ! Group 1
354:            J1=ATOMLISTS(T1,1,J2)  ! Actual atom index247:            J1=ATOMLISTS(T1,1,J2)  ! Actual atom index
355:            J13=3*(J1-1)248:            J13=3*(J1-1)
356:            DO I=1,3249:            DO I=1,3
357:               GRAD(J13+I)=0.0D0250:               V(J13+I)=V(J13+I)+VTEMP(J1,I)
 251:               IF(STRESST) THEN ! Accumulate overall stress tensor
 252:                  DO J=I,3
 253:                     STRESS(0,I,J) = STRESS(0,I,J) + STRESS(J1,I,J)
 254:                     STRESS(J1,J,I) = STRESS(J1,I,J) ! Impose symmetry
 255:                  ENDDO
 256:               ENDIF
358:            ENDDO257:            ENDDO
359:         ENDDO258:         ENDDO
360:      ENDDO259:      ENDDO
361:      !260:      !
362:      IF(STRESST) THEN261:      IF(STRESST) THEN
363:         DO I=1,3262:         DO I=1,3
364:            DO J=I,3263:            DO J=I,3
 264:               STRESS(0,I,J) = 0.5D0*STRESS(0,I,J) ! We double-counted!
365:               STRESS(0,J,I) = STRESS(0,I,J) ! Impose symmetry265:               STRESS(0,J,I) = STRESS(0,I,J) ! Impose symmetry
366:            ENDDO266:            ENDDO
367:         ENDDO267:         ENDDO
368:      ENDIF268:      ENDIF
369:      !269:      !
370:      ! IF(HESST) THEN270:      ! TEST PERATOM ENERGIES AND GRADIENTS
371:      !    write(*,*) 'Hessian 6x6 block:'271:      ! DO T1=1,NSPECIES(0)
372:      !    DO J1=1,6272:      !    DO G1=1,2
373:      !       WRITE(*,'(6(1X,F10.5))') (HESS(J1,J2), J2=1,6)273:      !       WRITE(*,'(A,3(1X,I5))') "BGUPTA> T,G,NAT(T,G)=", &
374:      !    ENDDO              274:      !            T1, G1, ATOMLISTS(T1,G1,0)
375:      ! ENDIF275:      !       DO J2=1,ATOMLISTS(T1,G1,0)
 276:      !          J1=ATOMLISTS(T1,G1,J2)
 277:      !          J13=3*(J1-1)
 278:      !          WRITE(*,'(A,1X, I4, 4(1X,F12.6))') "BGUPTA>", &
 279:      !               J1, VT(J1), V(J13+1), V(J13+2), V(J13+3)
 280:      !       ENDDO
 281:      !    ENDDO
 282:      ! ENDDO
376:      !283:      !
377:   ENDIF 284:   ENDIF
378:   !285:   !
379:   RETURN286:   RETURN
380:   !287:   !
381: END SUBROUTINE MGUPTA288: END SUBROUTINE MGUPTA


r29119/parse_pot_params.f90 2015-11-17 23:33:06.740595148 +0000 r29118/parse_pot_params.f90 2015-11-17 23:33:08.112613547 +0000
148:     !148:     !
149:     RETURN149:     RETURN
150:     !150:     !
151:   END SUBROUTINE PARSE_MLJ_PARAMS151:   END SUBROUTINE PARSE_MLJ_PARAMS
152:   !152:   !
153:   SUBROUTINE PARSE_MGUPTA_PARAMS(DATA_UNIT)153:   SUBROUTINE PARSE_MGUPTA_PARAMS(DATA_UNIT)
154:     !154:     !
155:     ! Parse parameters for Multicomponent Gupta system.155:     ! Parse parameters for Multicomponent Gupta system.
156:     !156:     !
157:     USE POT_PARAMS, ONLY : MGUPTA_A, MGUPTA_XI, MGUPTA_P, &157:     USE POT_PARAMS, ONLY : MGUPTA_A, MGUPTA_XI, MGUPTA_P, &
158:          MGUPTA_Q, MGUPTA_R0, MGUPTA_M2Q, MGUPTA_XI_SQ, &158:          MGUPTA_Q, MGUPTA_R0
159:          MGUPTA_MP_DIVBY_R0, MGUPTA_M2Q_DIVBY_R0 
160:     !159:     !
161:     IMPLICIT NONE160:     IMPLICIT NONE
162:     !161:     !
163:     INTEGER, INTENT(IN) :: DATA_UNIT162:     INTEGER, INTENT(IN) :: DATA_UNIT
164:     !163:     !
165:     INTEGER :: N0,N1,N2,N3,NONTYPEA,NTYPEA164:     INTEGER :: N0,N1,N2,N3,NONTYPEA,NTYPEA
166:     !165:     !
167:     ! First check the number of species from the first argument166:     ! First check the number of species from the first argument
168:     CALL READI(N0)167:     CALL READI(N0)
169:     IF (NSPECIES(0) /= N0) THEN168:     IF (NSPECIES(0) /= N0) THEN
171:             'parse_pot_params> Inconsistent species count for MGUPTA!'170:             'parse_pot_params> Inconsistent species count for MGUPTA!'
172:        STOP171:        STOP
173:     ENDIF172:     ENDIF
174:     !173:     !
175:     IF(ALLOCATED(MGUPTA_A)) DEALLOCATE(MGUPTA_A)174:     IF(ALLOCATED(MGUPTA_A)) DEALLOCATE(MGUPTA_A)
176:     IF(ALLOCATED(MGUPTA_P)) DEALLOCATE(MGUPTA_P)175:     IF(ALLOCATED(MGUPTA_P)) DEALLOCATE(MGUPTA_P)
177:     IF(ALLOCATED(MGUPTA_Q)) DEALLOCATE(MGUPTA_Q)176:     IF(ALLOCATED(MGUPTA_Q)) DEALLOCATE(MGUPTA_Q)
178:     IF(ALLOCATED(MGUPTA_XI)) DEALLOCATE(MGUPTA_XI)177:     IF(ALLOCATED(MGUPTA_XI)) DEALLOCATE(MGUPTA_XI)
179:     IF(ALLOCATED(MGUPTA_R0)) DEALLOCATE(MGUPTA_R0)178:     IF(ALLOCATED(MGUPTA_R0)) DEALLOCATE(MGUPTA_R0)
180:     !179:     !
181:     IF(ALLOCATED(MGUPTA_M2Q)) DEALLOCATE(MGUPTA_M2Q) 
182:     IF(ALLOCATED(MGUPTA_XI_SQ)) DEALLOCATE(MGUPTA_XI_SQ) 
183:     IF(ALLOCATED(MGUPTA_MP_DIVBY_R0)) DEALLOCATE(MGUPTA_MP_DIVBY_R0) 
184:     IF(ALLOCATED(MGUPTA_M2Q_DIVBY_R0)) & 
185:          DEALLOCATE(MGUPTA_M2Q_DIVBY_R0) 
186:     ! 
187:     ALLOCATE(MGUPTA_A(N0,N0),MGUPTA_XI(N0,N0))180:     ALLOCATE(MGUPTA_A(N0,N0),MGUPTA_XI(N0,N0))
188:     ALLOCATE(MGUPTA_P(N0,N0),MGUPTA_Q(N0,N0),MGUPTA_R0(N0,N0))181:     ALLOCATE(MGUPTA_P(N0,N0),MGUPTA_Q(N0,N0),MGUPTA_R0(N0,N0))
189:     ALLOCATE(MGUPTA_M2Q(N0,N0),MGUPTA_XI_SQ(N0,N0)) 
190:     ALLOCATE(MGUPTA_MP_DIVBY_R0(N0,N0)) 
191:     ALLOCATE(MGUPTA_M2Q_DIVBY_R0(N0,N0)) 
192:     !182:     !
193:     CALL READF(MGUPTA_A(1,1))183:     CALL READF(MGUPTA_A(1,1))
194:     CALL READF(MGUPTA_P(1,1))184:     CALL READF(MGUPTA_P(1,1))
195:     CALL READF(MGUPTA_Q(1,1))185:     CALL READF(MGUPTA_Q(1,1))
196:     CALL READF(MGUPTA_XI(1,1))186:     CALL READF(MGUPTA_XI(1,1))
197:     CALL READF(MGUPTA_R0(1,1))187:     CALL READF(MGUPTA_R0(1,1))
198:     !188:     !
199:     NONTYPEA=0189:     NONTYPEA=0
200:     DO N1 = 1,NSPECIES(0) 190:     DO N1 = 1,NSPECIES(0) 
201:        !191:        !
246:           ELSE236:           ELSE
247:              WRITE(MYUNIT,'(A)') &237:              WRITE(MYUNIT,'(A)') &
248:                   'parse_MGupta_params> Missing ''MGUPTA'' header!'238:                   'parse_MGupta_params> Missing ''MGUPTA'' header!'
249:              STOP239:              STOP
250:           ENDIF240:           ENDIF
251:           !241:           !
252:        ENDDO242:        ENDDO
253:        !243:        !
254:     ENDDO244:     ENDDO
255:     !245:     !
256:     ! Also pre-compute some auxilliary constants for speed 
257:     DO N1 = 1,NSPECIES(0)  
258:        DO N2 = N1,NSPECIES(0) 
259:           ! 
260:           MGUPTA_XI_SQ(N1,N2) = MGUPTA_XI(N1,N2)**2 
261:           MGUPTA_XI_SQ(N2,N1) = MGUPTA_XI_SQ(N1,N2) 
262:           ! 
263:           MGUPTA_MP_DIVBY_R0(N1,N2) = & 
264:                -MGUPTA_P(N1,N2)/MGUPTA_R0(N1,N2) 
265:           MGUPTA_MP_DIVBY_R0(N2,N1) = MGUPTA_MP_DIVBY_R0(N1,N2) 
266:           ! 
267:           MGUPTA_M2Q(N1,N2) = -2.0D0*MGUPTA_Q(N1,N2) 
268:           MGUPTA_M2Q(N2,N1) = MGUPTA_M2Q(N1,N2) 
269:           ! 
270:           MGUPTA_M2Q_DIVBY_R0(N1,N2) = & 
271:                MGUPTA_M2Q(N1,N2)/MGUPTA_R0(N1,N2) 
272:           MGUPTA_M2Q_DIVBY_R0(N2,N1) = MGUPTA_M2Q_DIVBY_R0(N1,N2) 
273:           ! 
274:        ENDDO 
275:     ENDDO 
276:     ! 
277:     ! Atom number for species 1 has not yet been specified.        246:     ! Atom number for species 1 has not yet been specified.        
278:     NTYPEA = NATOMS - NONTYPEA247:     NTYPEA = NATOMS - NONTYPEA
279:     NSPECIES(1) = NTYPEA248:     NSPECIES(1) = NTYPEA
280:     !                                                              249:     !                                                              
281:     WRITE(MYUNIT,'(A,I4,A)') &250:     WRITE(MYUNIT,'(A,I4,A)') &
282:          'parse_MGupta_params> Gupta system with', &251:          'parse_MGupta_params> Gupta system with', &
283:          NSPECIES(0),' species.' 252:          NSPECIES(0),' species.' 
284:     WRITE(MYUNIT, '(A)') &253:     WRITE(MYUNIT, '(A)') &
285:          'parse_MGupta_params> Atom count for each species:'254:          'parse_MGupta_params> Atom count for each species:'
286:     DO N1=1,NSPECIES(0)255:     DO N1=1,NSPECIES(0)


r29119/potential.f90 2015-11-17 23:33:07.116600189 +0000 r29118/potential.f90 2015-11-17 23:33:08.488618591 +0000
 70:    REAL(REAL64)               :: PTGPSAVE(3, 3, 2*CSMGPINDEX), CSMNORMSAVE, ENERGY, VNEW(3*NATOMS) 70:    REAL(REAL64)               :: PTGPSAVE(3, 3, 2*CSMGPINDEX), CSMNORMSAVE, ENERGY, VNEW(3*NATOMS)
 71:    INTEGER(INT32)             :: BRUN 71:    INTEGER(INT32)             :: BRUN
 72:    INTEGER(INT32)             :: NQTOT 72:    INTEGER(INT32)             :: NQTOT
 73:    LOGICAL                    :: SOCOUPLE, GUIDECHANGET, CSMDOGUIDET, COMTEST 73:    LOGICAL                    :: SOCOUPLE, GUIDECHANGET, CSMDOGUIDET, COMTEST
 74: ! khs26 > AMBER12 energy decomposition, type defined in AMBER12 interface 74: ! khs26 > AMBER12 energy decomposition, type defined in AMBER12 interface
 75:    TYPE(POT_ENE_REC_C)        :: AMBER12_ENERGY_DECOMP 75:    TYPE(POT_ENE_REC_C)        :: AMBER12_ENERGY_DECOMP
 76: ! Used only in commented sections 76: ! Used only in commented sections
 77: !   REAL(REAL64)               :: XG, YG, ZG, GRADLJ(3*NATOMS), EREALLJ, GRADMF(3*NATOMS), EREALMF, TERMLJ, TERMMF & 77: !   REAL(REAL64)               :: XG, YG, ZG, GRADLJ(3*NATOMS), EREALLJ, GRADMF(3*NATOMS), EREALMF, TERMLJ, TERMMF &
 78: !                                 GRADDUM(3*NATOMS), GRADDUM1(3*NATOMS), GRADDUM2(3*NATOMS), EPLUS, EMINUS, DIFF 78: !                                 GRADDUM(3*NATOMS), GRADDUM1(3*NATOMS), GRADDUM2(3*NATOMS), EPLUS, EMINUS, DIFF
 79: !                                 EDUM1, EDUM2 79: !                                 EDUM1, EDUM2
 80: !DS656> GRADIENT/HESSIAN TESTING 
 81:    DOUBLE PRECISION :: GPLUS(3*NATOMS), GMINUS(3*NATOMS), DIFF, & 
 82:         EPLUS, EMINUS 
 83:     
 84: ! TODO: Delete common blocks 80: ! TODO: Delete common blocks
 85:    COMMON /CO/ COMPON 81:    COMMON /CO/ COMPON
 86:    COMMON /FAIL/ FTEST 82:    COMMON /FAIL/ FTEST
 87:    COMMON /EV/ EVAP, EVAPREJECT 83:    COMMON /EV/ EVAP, EVAPREJECT
 88:    COMMON /GD/ GUIDECHANGET, GUIDET, CSMDOGUIDET 84:    COMMON /GD/ GUIDECHANGET, GUIDET, CSMDOGUIDET
 89:    COMMON /CSMAVVAL/ AVVAL, CSMRMS, CSMIT 85:    COMMON /CSMAVVAL/ AVVAL, CSMRMS, CSMIT
 90:    COMMON /TOT/ NQTOT 86:    COMMON /TOT/ NQTOT
 91:  87: 
 92:    GUIDECHANGET=.FALSE. 88:    GUIDECHANGET=.FALSE.
 93: ! 89: !
598:    ELSE IF (FST) THEN594:    ELSE IF (FST) THEN
599:       CALL RAD(X, GRAD, EREAL, GRADT)595:       CALL RAD(X, GRAD, EREAL, GRADT)
600:       CALL FS(X, GRAD, EREAL, GRADT)596:       CALL FS(X, GRAD, EREAL, GRADT)
601: 597: 
602:    ELSE IF (BGUPTAT) THEN598:    ELSE IF (BGUPTAT) THEN
603:       CALL RAD(X, GRAD, EREAL, GRADT)599:       CALL RAD(X, GRAD, EREAL, GRADT)
604:       CALL BGUPTA(X, GRAD, EREAL, GRADT)600:       CALL BGUPTA(X, GRAD, EREAL, GRADT)
605: 601: 
606:    ELSE IF (MGUPTAT) THEN602:    ELSE IF (MGUPTAT) THEN
607:       CALL RAD(X, GRAD, EREAL, GRADT)603:       CALL RAD(X, GRAD, EREAL, GRADT)
608:       CALL MGUPTA(X, GRAD, EREAL, GRADT, SECT, .FALSE.)604:       CALL MGUPTA(X, GRAD, EREAL, GRADT, .FALSE.)
609:       ! Testing...605: 
610:       !ds656> Test gradient and Hessian 
611:       IF (SECT.AND..FALSE.) THEN 
612:          DIFF=1.0D-5 
613:          PRINT*,SECT 
614:          PRINT*,'analytic and numerical gradients:' 
615:          DO J1=1,3*NATOMS 
616:             X(J1)=X(J1)+DIFF 
617:             CALL MGUPTA(X,GPLUS,EPLUS,.FALSE.,.FALSE.,.false.) 
618:             X(J1)=X(J1)-2.0D0*DIFF 
619:             CALL MGUPTA(X,GMINUS,EMINUS,.FALSE.,.FALSE.,.false.) 
620:             X(J1)=X(J1)+DIFF 
621:             IF ((ABS(GRAD(J1)).NE.0.0D0).AND.(100.0D0*(GRAD(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GRAD(J1).GT.1.0D0)) THEN 
622:                WRITE(*,'(I5,2F20.10)') J1,GRAD(J1),(EPLUS-EMINUS)/(2.0D0*DIFF) 
623:             ENDIF 
624:          ENDDO 
625:          PRINT*,'analytic and numerical second derivatives:' 
626:          DO J1=1,3*NATOMS 
627:             X(J1)=X(J1)+DIFF 
628:             CALL MGUPTA(X,GPLUS,EPLUS,.TRUE.,.FALSE.,.false.) 
629:             X(J1)=X(J1)-2.0D0*DIFF 
630:             CALL MGUPTA(X,GMINUS,EMINUS,.TRUE.,.FALSE.,.false.) 
631:             X(J1)=X(J1)+DIFF 
632:             DO J2=1,3*NATOMS 
633:                IF ((ABS(HESS(J1,J2)).NE.0.0D0).AND. & 
634:                     (ABS(100.0D0*(HESS(J1,J2)-(GPLUS(J2)-GMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D0)) THEN 
635:                   WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(GPLUS(J2)-GMINUS(J2))/(2.0D0*DIFF),'   X' 
636:                ELSE 
637:                   WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(GPLUS(J2)-GMINUS(J2))/(2.0D0*DIFF) 
638:                ENDIF 
639:             ENDDO 
640:          ENDDO 
641:          STOP 
642:       ENDIF 
643:       ! <ds656 ...end of testing. 
644:    ELSE IF (GLJY) THEN606:    ELSE IF (GLJY) THEN
645:       CALL RAD(X, GRAD, EREAL, GRADT)607:       CALL RAD(X, GRAD, EREAL, GRADT)
646:       CALL GLJYPOT(X, GRAD, EREAL, GRADT)608:       CALL GLJYPOT(X, GRAD, EREAL, GRADT)
647: 609: 
648:    ELSE IF (GUPTAT) THEN610:    ELSE IF (GUPTAT) THEN
649:       CALL RAD(X, GRAD, EREAL, GRADT)611:       CALL RAD(X, GRAD, EREAL, GRADT)
650:       CALL GUPTA(X, GRAD, EREAL, GRADT)612:       CALL GUPTA(X, GRAD, EREAL, GRADT)
651: 613: 
652:    ELSE IF (NATBT) THEN614:    ELSE IF (NATBT) THEN
653:       CALL RAD(X, GRAD, EREAL, GRADT)615:       CALL RAD(X, GRAD, EREAL, GRADT)
1104:       EREAL=EREAL - PFORCE*(X(3*(PATOM1-1)+3)-X(3*(PATOM2-1)+3))1066:       EREAL=EREAL - PFORCE*(X(3*(PATOM1-1)+3)-X(3*(PATOM2-1)+3))
1105:       GRAD(3*(PATOM1-1)+3)=GRAD(3*(PATOM1-1)+3) - PFORCE1067:       GRAD(3*(PATOM1-1)+3)=GRAD(3*(PATOM1-1)+3) - PFORCE
1106:       GRAD(3*(PATOM2-1)+3)=GRAD(3*(PATOM2-1)+3) + PFORCE1068:       GRAD(3*(PATOM2-1)+3)=GRAD(3*(PATOM2-1)+3) + PFORCE
1107:    END IF1069:    END IF
1108: 1070: 
1109:    IF (GCBHT) THEN1071:    IF (GCBHT) THEN
1110: ! Need a negative sign because the exponential is -(V-N mu)/kT1072: ! Need a negative sign because the exponential is -(V-N mu)/kT
1111:       EREAL=EREAL - NATOMS*GCMU1073:       EREAL=EREAL - NATOMS*GCMU
1112:    END IF1074:    END IF
1113: 1075: 
 1076: ! ds656> If relative chemical potentials are present, then we tack
 1077: !        them on to the energy for semi-grand canonical BH.
 1078:    IF (SEMIGRAND_MUT) THEN      
 1079:       DO J1=2, NSPECIES(0)
 1080:          EREAL = EREAL - NSPECIES(J1)*SEMIGRAND_MU(J1)
 1081:       ENDDO
 1082:    ENDIF
 1083:    
1114:    IF (FIELDT) THEN1084:    IF (FIELDT) THEN
1115:       IF (CENT) THEN1085:       IF (CENT) THEN
1116:          CALL FDM(X, GRAD, EREAL, GRADT)1086:          CALL FDM(X, GRAD, EREAL, GRADT)
1117:       ELSE1087:       ELSE
1118:          CALL FD(X, GRAD, EREAL, GRADT)1088:          CALL FD(X, GRAD, EREAL, GRADT)
1119:       END IF1089:       END IF
1120:    ELSE IF (CPMD .AND. (NPCALL == 1)) THEN1090:    ELSE IF (CPMD .AND. (NPCALL == 1)) THEN
1121:       CALL SYSTEM(' sed -e "s/DUMMY/RESTART WAVEFUNCTION GEOFILE LATEST/" ' //  SYS(1:LSYS) // ' > temp ')1091:       CALL SYSTEM(' sed -e "s/DUMMY/RESTART WAVEFUNCTION GEOFILE LATEST/" ' //  SYS(1:LSYS) // ' > temp ')
1122:       CALL SYSTEM(' mv temp ' // SYS(1:LSYS) // '.restart')1092:       CALL SYSTEM(' mv temp ' // SYS(1:LSYS) // '.restart')
1123:    ELSE IF (CPMD .AND. (.NOT. SCT)) THEN1093:    ELSE IF (CPMD .AND. (.NOT. SCT)) THEN


r29119/pot_params.f90 2015-11-17 23:33:06.924597613 +0000 r29118/pot_params.f90 2015-11-17 23:33:08.300616068 +0000
 26: MODULE POT_PARAMS 26: MODULE POT_PARAMS
 27:   ! 27:   !
 28:   IMPLICIT NONE 28:   IMPLICIT NONE
 29:   ! 29:   !
 30:   ! --- Multicomponent Lennard-Jones parameters ---------------- 30:   ! --- Multicomponent Lennard-Jones parameters ----------------
 31:   DOUBLE PRECISION, ALLOCATABLE :: MLJ_EPS(:,:), MLJ_SIG(:,:)  31:   DOUBLE PRECISION, ALLOCATABLE :: MLJ_EPS(:,:), MLJ_SIG(:,:) 
 32:   ! ------------------------------------------------------------ 32:   ! ------------------------------------------------------------
 33:   ! 33:   !
 34:   ! --- Multicomponent Gupta parameters ------------------------ 34:   ! --- Multicomponent Gupta parameters ------------------------
 35:   DOUBLE PRECISION, ALLOCATABLE :: MGUPTA_A(:,:), MGUPTA_XI(:,:), & 35:   DOUBLE PRECISION, ALLOCATABLE :: MGUPTA_A(:,:), MGUPTA_XI(:,:), &
 36:        MGUPTA_P(:,:), MGUPTA_Q(:,:), MGUPTA_R0(:,:), & 36:        MGUPTA_P(:,:), MGUPTA_Q(:,:), MGUPTA_R0(:,:)
 37:        MGUPTA_M2Q(:,:), MGUPTA_XI_SQ(:,:), & 
 38:        MGUPTA_MP_DIVBY_R0(:,:), MGUPTA_M2Q_DIVBY_R0(:,:) 
 39:   ! ------------------------------------------------------------ 37:   ! ------------------------------------------------------------
 40:   ! 38:   !
 41:   ! --- Multicomponent Sutton-Chen parameters ------------------ 39:   ! --- Multicomponent Sutton-Chen parameters ------------------
 42:   INTEGER, ALLOCATABLE :: MSC_N(:,:), MSC_M(:,:) 40:   INTEGER, ALLOCATABLE :: MSC_N(:,:), MSC_M(:,:)
 43:   DOUBLE PRECISION, ALLOCATABLE :: MSC_EPS(:,:), MSC_A(:,:), & 41:   DOUBLE PRECISION, ALLOCATABLE :: MSC_EPS(:,:), MSC_A(:,:), &
 44:        MSC_C(:) 42:        MSC_C(:)
 45:   DOUBLE PRECISION, ALLOCATABLE :: CUTA_REP(:,:),CUTB_REP(:,:), & 43:   DOUBLE PRECISION, ALLOCATABLE :: CUTA_REP(:,:),CUTB_REP(:,:), &
 46:        CUTC_REP(:,:) 44:        CUTC_REP(:,:)
 47:   DOUBLE PRECISION, ALLOCATABLE :: CUTA_ATT(:,:),CUTB_ATT(:,:), & 45:   DOUBLE PRECISION, ALLOCATABLE :: CUTA_ATT(:,:),CUTB_ATT(:,:), &
 48:        CUTC_ATT(:,:) 46:        CUTC_ATT(:,:)


r29119/quench.F 2015-11-17 23:33:07.340603195 +0000 r29118/quench.F 2015-11-17 23:33:08.704621486 +0000
319:             CALL CPU_TIME(T_LBFGS_START)319:             CALL CPU_TIME(T_LBFGS_START)
320:             CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)320:             CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)
321:             CALL CPU_TIME(T_LBFGS_END)321:             CALL CPU_TIME(T_LBFGS_END)
322:             IF (FEBHT) THEN322:             IF (FEBHT) THEN
323:                WRITE(MYUNIT, '(A, F10.3)') 'Time to minimise (s):', T_LBFGS_END - T_LBFGS_START323:                WRITE(MYUNIT, '(A, F10.3)') 'Time to minimise (s):', T_LBFGS_END - T_LBFGS_START
324:             ENDIF324:             ENDIF
325:          ENDIF325:          ENDIF
326:          IF (EVAPREJECT) RETURN326:          IF (EVAPREJECT) RETURN
327:       ENDIF327:       ENDIF
328: 328: 
329: ! ds656> If (relative) chemical potentials are present, then we tack 
330: !        them on to the energy for semi-grand canonical BH. 
331:       IF (SEMIGRAND_MUT) THEN      
332:          DO J1=2, NSPECIES(0) 
333:             EREAL = EREAL - NSPECIES(J1)*SEMIGRAND_MU(J1) 
334:          ENDDO 
335:       ENDIF 
336:  
337:       IF (FEBHT .AND. CFLAG) THEN329:       IF (FEBHT .AND. CFLAG) THEN
338:       ! Calculate the free energy330:       ! Calculate the free energy
339:          NUM_ZERO_EVS=6331:          NUM_ZERO_EVS=6
340:          IF (NATOMS.EQ.2) NUM_ZERO_EVS=5332:          IF (NATOMS.EQ.2) NUM_ZERO_EVS=5
341:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)333:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
342:          ALLOCATE(HESS(3*NATOMS,3*NATOMS))334:          ALLOCATE(HESS(3*NATOMS,3*NATOMS))
343:          CALL POTENTIAL(P,GRAD,EREAL,.TRUE.,.TRUE.)335:          CALL POTENTIAL(P,GRAD,EREAL,.TRUE.,.TRUE.)
344:          CALL MASSWT()336:          CALL MASSWT()
345: #ifdef __SPARSE337: #ifdef __SPARSE
346:          IF (SPARSET) THEN338:          IF (SPARSET) THEN


r29119/stress.f90 2015-11-17 23:33:07.532605768 +0000 r29118/stress.f90 2015-11-17 23:33:08.896624063 +0000
 13:   ! 13:   !
 14:   IF(.NOT.STRESST) THEN 14:   IF(.NOT.STRESST) THEN
 15:      WRITE(MYUNIT,'(A)') & 15:      WRITE(MYUNIT,'(A)') &
 16:           'calc_stress> Should not be here!' 16:           'calc_stress> Should not be here!'
 17:      RETURN 17:      RETURN
 18:   ELSE 18:   ELSE
 19:      STRESS(:,:,:) = 0.0D0 ! Initialise 19:      STRESS(:,:,:) = 0.0D0 ! Initialise
 20:   ENDIF 20:   ENDIF
 21:   ! 21:   !
 22:   IF(MGUPTAT) THEN 22:   IF(MGUPTAT) THEN
 23:      CALL MGUPTA(X, GRAD, EPOT, .TRUE., .FALSE., .TRUE.) 23:      CALL MGUPTA(X, GRAD, EPOT, .TRUE., .TRUE.)
 24:   ELSE 24:   ELSE
 25:      WRITE(MYUNIT,'(A)') & 25:      WRITE(MYUNIT,'(A)') &
 26:           'calc_stress> Stress calculation not implemented for current potential.' 26:           'calc_stress> Cannot compute stress for current pot.'
 27:   ENDIF 27:   ENDIF
 28:   ! 28:   !
 29:   IF(MIEFT) CALL MIEF(X,GRAD,EPOT,.TRUE.,.TRUE.) 29:   IF(MIEFT) CALL MIEF(X,GRAD,EPOT,.TRUE.,.TRUE.)
 30:   ! 30:   !
 31:   WRITE(MYUNIT,'(A)') 'stress> Overall stress tensor:' 31:   WRITE(MYUNIT,'(A)') 'stress> Overall stress tensor:'
 32:   DO I=1,3 32:   DO I=1,3
 33:      WRITE(MYUNIT,'(3(1X,E15.8))') (STRESS(0,I,J),J=1,3) 33:      WRITE(MYUNIT,'(3(1X,E15.8))') (STRESS(0,I,J),J=1,3)
 34:   ENDDO 34:   ENDDO
 35:   ! 35:   !
 36:   RETURN 36:   RETURN


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0