hdiff output

r31119/atom_label_swaps.f90 2016-09-06 15:30:37.470425233 +0100 r31118/atom_label_swaps.f90 2016-09-06 15:30:40.694468265 +0100
 52:      DO J=1,NSPECIES(0) 52:      DO J=1,NSPECIES(0)
 53:         IF(J==I2) CYCLE 53:         IF(J==I2) CYCLE
 54:         DO K=1,ATOMLISTS(J,1,0) ! only group 1 54:         DO K=1,ATOMLISTS(J,1,0) ! only group 1
 55:            LIST(0) = LIST(0) + 1 55:            LIST(0) = LIST(0) + 1
 56:            LIST(LIST(0)) = ATOMLISTS(J,1,K)  56:            LIST(LIST(0)) = ATOMLISTS(J,1,K) 
 57:         ENDDO 57:         ENDDO
 58:      ENDDO 58:      ENDDO
 59:      I2 = INT(DPRAND()*DBLE(LIST(0))) + 1 ! Forgot +1 59:      I2 = INT(DPRAND()*DBLE(LIST(0))) + 1 ! Forgot +1
 60:      I2 = LIST(I2) 60:      I2 = LIST(I2)
 61:      ! 61:      !
 62:      CALL SWAP_LABELS(I1,I2,NP) 62:      CALL SWAP_LABELS(I1,I2)
 63:      ! 63:      !
 64:      NQTOT = NQTOT + 1 64:      NQTOT = NQTOT + 1
 65:      NQ(NP) = NQ(NP) + 1 65:      NQ(NP) = NQ(NP) + 1
 66:      CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC) 66:      CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)
 67:      IF (MOD(N-1,PRTFRQ).EQ.0) THEN 67:      IF (MOD(N-1,PRTFRQ).EQ.0) THEN
 68:         IF (NPAR.GT.1) THEN 68:         IF (NPAR.GT.1) THEN
 69:            WRITE(MYUNIT,'(A,I1,A,I10,A,F20.10,A,I5,A,G12.5,A,F11.1)') & 69:            WRITE(MYUNIT,'(A,I1,A,I10,A,F20.10,A,I5,A,G12.5,A,F11.1)') &
 70:                 '[',NP,']Qu ',NQ(NP),' E=',POTEL,' steps=',ITER, & 70:                 '[',NP,']Qu ',NQ(NP),' E=',POTEL,' steps=',ITER, &
 71:                 ' RMS=',RMS,' t=',TIME 71:                 ' RMS=',RMS,' t=',TIME
 72:         ELSE 72:         ELSE
 84:            XMIN(1:3*NATOMS)=COORDS(1:3*NATOMS, NP) 84:            XMIN(1:3*NATOMS)=COORDS(1:3*NATOMS, NP)
 85:            ATOMLISTS_MIN(:,:,:) = ATOMLISTS(:,:,:) 85:            ATOMLISTS_MIN(:,:,:) = ATOMLISTS(:,:,:)
 86:            INVATOMLISTS_MIN(:,:) = INVATOMLISTS(:,:) 86:            INVATOMLISTS_MIN(:,:) = INVATOMLISTS(:,:)
 87:            WRITE(MYUNIT,'(A,F20.10,A,I8,A,F10.8)') & 87:            WRITE(MYUNIT,'(A,F20.10,A,I8,A,F10.8)') &
 88:                 'uswaps> New EMIN= ',EMIN,' on swap', N,& 88:                 'uswaps> New EMIN= ',EMIN,' on swap', N,&
 89:                 ' at T=', TEMP 89:                 ' at T=', TEMP
 90:            IF(HIT) RETURN 90:            IF(HIT) RETURN
 91:         ENDIF 91:         ENDIF
 92:      ELSE ! Reject. Revert labels and coordinates. 92:      ELSE ! Reject. Revert labels and coordinates.
 93:         NREJ = NREJ + 1 93:         NREJ = NREJ + 1
 94:         CALL SWAP_LABELS(I1,I2,NP) 94:         CALL SWAP_LABELS(I1,I2)
 95:         COORDS(1:3*NATOMS,NP) = X0(1:3*NATOMS) 95:         COORDS(1:3*NATOMS,NP) = X0(1:3*NATOMS)
 96:         POTEL = E0 96:         POTEL = E0
 97:      ENDIF 97:      ENDIF
 98:      ! 98:      !
 99:      IF(MOD(N,LSWAPS_NUP)==0) THEN 99:      IF(MOD(N,LSWAPS_NUP)==0) THEN
100:         IF(LSWAPS_TFACTOR <= 1.0D0) THEN 100:         IF(LSWAPS_TFACTOR <= 1.0D0) THEN 
101:            ! Simulated annealing: always decrease temperature. 101:            ! Simulated annealing: always decrease temperature. 
102:            TEMP = TEMP*LSWAPS_TFACTOR102:            TEMP = TEMP*LSWAPS_TFACTOR
103:         ELSE ! Strive for ACC/REJ ratio of 0.5.103:         ELSE ! Strive for ACC/REJ ratio of 0.5.
104:            ! This will have to be recoded for desired ACC/REJ 104:            ! This will have to be recoded for desired ACC/REJ 


r31119/atomlists.f90 2016-09-06 15:30:37.730428697 +0100 r31118/atomlists.f90 2016-09-06 15:30:40.942471574 +0100
 67:      DO I=1,NATOMS 67:      DO I=1,NATOMS
 68:         ATMASS(I) = SPECMASS(INVATOMLISTS(I,1)) 68:         ATMASS(I) = SPECMASS(INVATOMLISTS(I,1))
 69:         WRITE(MYUNIT,*) ATMASS(I) 69:         WRITE(MYUNIT,*) ATMASS(I)
 70:      ENDDO 70:      ENDDO
 71:   ENDIF 71:   ENDIF
 72:   ! 72:   !
 73:   RETURN 73:   RETURN
 74:   ! 74:   !
 75: END SUBROUTINE RESET_ATOMLISTS 75: END SUBROUTINE RESET_ATOMLISTS
 76: ! 76: !
 77: SUBROUTINE RESET_LABELS(IP) 
 78:   ! 
 79:   USE COMMONS, ONLY : LABELS, NATOMS, NATOMSALLOC, NSPECIES, & 
 80:        SPECMASST, SPECMASS, ATMASS, MYUNIT 
 81:   ! 
 82:   IMPLICIT NONE 
 83:   ! 
 84:   INTEGER, INTENT(IN) :: IP 
 85:   INTEGER :: I,J,K 
 86:   ! 
 87:   LABELS(:,IP) = 0 
 88:   ! 
 89:   ! Separate atoms by type/species and make them all of group IGROUP 
 90:   K = 0 
 91:   DO I = 1, NSPECIES(0) ! span all species 
 92:      DO J = 1, NSPECIES(I) ! span all atoms of each species 
 93:         K = K+1 ! increment atom index 
 94:         LABELS(K, IP) = I 
 95:      ENDDO 
 96:   ENDDO 
 97:   ! 
 98:   ! Sanity check 
 99:   IF(NATOMS /= K) THEN 
100:      WRITE(MYUNIT, '(A)') 'reset_labels> Inconsistent atom count!' 
101:      !WRITE(MYUNIT, *) NATOMS, K, NSPECIES(0) 
102:      STOP 
103:   ENDIF 
104:   ! 
105:   IF(SPECMASST) THEN 
106:      WRITE(MYUNIT,'(A)') & 
107:           'initialization> (Re)setting masses in accord with SPECMASS:' 
108:      DO I=1,NATOMS 
109:         ATMASS(I) = SPECMASS(LABELS(I,IP)) 
110:         WRITE(MYUNIT,*) ATMASS(I) 
111:      ENDDO 
112:   ENDIF   
113:   ! 
114:   RETURN 
115:   ! 
116: END SUBROUTINE RESET_LABELS 
117: ! 
118: ! ds656> Routine for resetting ATOMLISTS in accord with provided 77: ! ds656> Routine for resetting ATOMLISTS in accord with provided
119: ! list of atomic labels.  78: ! list of atomic labels. 
120: SUBROUTINE SET_ATOMLISTS(LABELS,IGROUP) 79: SUBROUTINE SET_ATOMLISTS(LABELS,IGROUP)
121:   ! 80:   !
122:   USE COMMONS, ONLY : ATOMLISTS, INVATOMLISTS, NATOMS, NSPECIES, & 81:   USE COMMONS, ONLY : ATOMLISTS, INVATOMLISTS, NATOMS, NSPECIES, &
123:        SPECMASST, SPECMASS, ATMASS, MYUNIT, LFLIPST, QALCSMODE 82:        SPECMASST, SPECMASS, ATMASS, MYUNIT, LFLIPST, QALCSMODE
124:   ! 83:   !
125:   IMPLICIT NONE 84:   IMPLICIT NONE
126:   ! 85:   !
127:   INTEGER, INTENT(IN) :: LABELS(NATOMS), IGROUP 86:   INTEGER, INTENT(IN) :: LABELS(NATOMS), IGROUP
156:         WRITE(MYUNIT,'(A, I2, A, I5, I5)') &115:         WRITE(MYUNIT,'(A, I2, A, I5, I5)') &
157:              'set_atomlists> Inconsistent counts for atom-type ', &116:              'set_atomlists> Inconsistent counts for atom-type ', &
158:              I,' ,these numbers should equal:', NSPECIES(I), &117:              I,' ,these numbers should equal:', NSPECIES(I), &
159:              TYPECOUNTS(I)118:              TYPECOUNTS(I)
160:         STOP119:         STOP
161:      ENDIF120:      ENDIF
162:   END DO121:   END DO
163:   !122:   !
164: END SUBROUTINE SET_ATOMLISTS123: END SUBROUTINE SET_ATOMLISTS
165: !124: !
166: SUBROUTINE SET_LABELS(L,IP) 
167:   ! 
168:   USE COMMONS, ONLY : LABELS, NATOMS, NSPECIES, & 
169:        SPECMASST, SPECMASS, ATMASS, MYUNIT, LFLIPST, QALCSMODE 
170:   ! 
171:   IMPLICIT NONE 
172:   ! 
173:   INTEGER, INTENT(IN) :: L(NATOMS), IP 
174:   INTEGER :: I,TYPECOUNTS(1:NSPECIES(0)) 
175:   ! 
176:   TYPECOUNTS(:) = 0 
177:   ! 
178:   DO I=1,NATOMS 
179:      LABELS(I,IP) = L(I) 
180:      TYPECOUNTS(L(I)) = TYPECOUNTS(L(I)) + 1 
181:      IF(SPECMASST) ATMASS(I) = SPECMASS(L(I)) 
182:   ENDDO 
183:   ! 
184:   DO I=1,NSPECIES(0) 
185:      IF(LFLIPST.OR.(QALCSMODE.GE.6)) THEN 
186:         NSPECIES(I) = TYPECOUNTS(I) 
187:      ELSE IF(NSPECIES(I) /= TYPECOUNTS(I)) THEN 
188:         WRITE(MYUNIT,'(A, I2, A, I5, I5)') & 
189:              'set_labels> Inconsistent counts for atom-type ', & 
190:              I,' ,these numbers should equal:', NSPECIES(I), & 
191:              TYPECOUNTS(I) 
192:         STOP 
193:      ENDIF 
194:   END DO 
195:   ! 
196:   RETURN 
197:   ! 
198: END SUBROUTINE SET_LABELS 
199: ! 
200: SUBROUTINE CHANGE_ATOM_GROUP(LIST,IGROUP)125: SUBROUTINE CHANGE_ATOM_GROUP(LIST,IGROUP)
201:   !126:   !
202:   ! Move all atoms in LIST from their current group to IGROUP 127:   ! Move all atoms in LIST from their current group to IGROUP 
203:   !128:   !
204:   USE COMMONS, ONLY : ATOMLISTS, INVATOMLISTS, NATOMS, MYUNIT129:   USE COMMONS, ONLY : ATOMLISTS, INVATOMLISTS, NATOMS, MYUNIT
205:   !130:   !
206:   IMPLICIT NONE131:   IMPLICIT NONE
207:   !132:   !
208:   INTEGER, INTENT(IN) :: LIST(0:NATOMS),IGROUP133:   INTEGER, INTENT(IN) :: LIST(0:NATOMS),IGROUP
209:   INTEGER :: I,IT,IG,IP,J,K134:   INTEGER :: I,IT,IG,IP,J,K
238:      ! Update inverse list of atom I163:      ! Update inverse list of atom I
239:      INVATOMLISTS(I,2) = IGROUP164:      INVATOMLISTS(I,2) = IGROUP
240:      INVATOMLISTS(I,3) = J165:      INVATOMLISTS(I,3) = J
241:      !166:      !
242:   ENDDO167:   ENDDO
243:   !168:   !
244:   RETURN169:   RETURN
245:   !170:   !
246: END SUBROUTINE CHANGE_ATOM_GROUP171: END SUBROUTINE CHANGE_ATOM_GROUP
247: !172: !
248: SUBROUTINE RESET_STOICHIOMETRY(IP)173: SUBROUTINE RESET_STOICHIOMETRY()
249:   !174:   !
250:   ! ds656> 31/07/2015175:   ! ds656> 31/07/2015
251:   ! Reset stoichiometry to the initial value, with the labels 176:   ! Reset stoichiometry to the initial value, with the labels 
252:   ! distributed randomly.177:   ! distributed randomly.
253:   !178:   !
254:   USE COMMONS, ONLY : NATOMS, MYUNIT, NSPECIES_INI, LABELS179:   USE COMMONS, ONLY : NATOMS, MYUNIT, NSPECIES_INI
255:   !180:   !
256:   IMPLICIT NONE181:   IMPLICIT NONE
257:   !182:   !
258:   INTEGER, INTENT(IN) :: IP183:   INTEGER :: I, J, K, LABELS(NATOMS)
259:   INTEGER :: I, J, K 
260:   !184:   !
261:   LABELS(:,IP) = 0185:   LABELS(:) = 0
262:   K = 0186:   K = 0
263:   DO I=1, NSPECIES_INI(0)187:   DO I=1, NSPECIES_INI(0)
264:      DO J= 1,NSPECIES_INI(I)188:      DO J= 1,NSPECIES_INI(I)
265:         K=K+1189:         K=K+1
266:         IF(K > NATOMS) THEN190:         IF(K > NATOMS) THEN
267:            WRITE(MYUNIT,'(A)') 'set_stoichiometry> WTF?!? Terminatig..'191:            WRITE(MYUNIT,'(A)') 'set_stoichiometry> WTF?!? Terminatig..'
268:            STOP192:            STOP
269:         ELSE193:         ELSE
270:            LABELS(K,IP) = I194:            LABELS(K) = I
271:         ENDIF195:         ENDIF
272:      ENDDO196:      ENDDO
273:   ENDDO197:   ENDDO
274:   !198:   !
275:   CALL KNUTH_SHUFFLE(NATOMS, LABELS(1:NATOMS,IP))199:   CALL KNUTH_SHUFFLE(NATOMS, LABELS)
276:   CALL SET_ATOMLISTS(LABELS(1:NATOMS,IP), 1) ! put all in group 1200:   CALL SET_ATOMLISTS(LABELS, 1) ! put all in group 1
277:   !201:   !
278:   RETURN202:   RETURN
279:   !203:   !
280: END SUBROUTINE RESET_STOICHIOMETRY204: END SUBROUTINE RESET_STOICHIOMETRY
281: !205: !
282: SUBROUTINE KNUTH_SHUFFLE(N,A)206: SUBROUTINE KNUTH_SHUFFLE(N,A)
283:   !207:   !
284:   ! ds656> Randomly permute integer vector A208:   ! ds656> Randomly permute integer vector A
285:   !209:   !
286:   INTEGER, INTENT(IN) :: N210:   INTEGER, INTENT(IN) :: N


r31119/commons.f90 2016-09-06 15:30:37.986432120 +0100 r31118/commons.f90 2016-09-06 15:30:41.250475687 +0100
400:       LOGICAL, ALLOCATABLE, DIMENSION(:) :: HARMONICFLIST  !  NATOMS400:       LOGICAL, ALLOCATABLE, DIMENSION(:) :: HARMONICFLIST  !  NATOMS
401:       INTEGER, ALLOCATABLE, DIMENSION(:) :: FROZENLIST  !  NATOMS401:       INTEGER, ALLOCATABLE, DIMENSION(:) :: FROZENLIST  !  NATOMS
402: 402: 
403:       INTEGER, ALLOCATABLE, DIMENSION(:) ::  IATNUM     403:       INTEGER, ALLOCATABLE, DIMENSION(:) ::  IATNUM     
404:       INTEGER,ALLOCATABLE :: ANV(:,:,:)404:       INTEGER,ALLOCATABLE :: ANV(:,:,:)
405:       DOUBLE PRECISION, ALLOCATABLE :: LJREP_BLN(:,:), LJATT_BLN(:,:), A_BLN(:), B_BLN(:), C_BLN(:), D_BLN(:), HYDRO_BLN(:)405:       DOUBLE PRECISION, ALLOCATABLE :: LJREP_BLN(:,:), LJATT_BLN(:,:), A_BLN(:), B_BLN(:), C_BLN(:), D_BLN(:), HYDRO_BLN(:)
406:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: TCOORDS406:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: TCOORDS
407:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: VT407:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: VT
408:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: VTVT408:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: VTVT
409:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: COORDS,COORDSO409:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: COORDS,COORDSO
410:       INTEGER, ALLOCATABLE, DIMENSION(:,:) :: LABELS, LABELSO ! <ds656 
411:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: FINISH410:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: FINISH
412:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) ::  VAT,  VATO    411:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) ::  VAT,  VATO    
413:       DOUBLE PRECISION, ALLOCATABLE :: SITE(:,:), RBUV(:,:), RBSTLA(:,:), STCHRG(:), DPMU(:)412:       DOUBLE PRECISION, ALLOCATABLE :: SITE(:,:), RBUV(:,:), RBSTLA(:,:), STCHRG(:), DPMU(:)
414:       DOUBLE PRECISION, ALLOCATABLE :: MSBCOORDS(:,:)413:       DOUBLE PRECISION, ALLOCATABLE :: MSBCOORDS(:,:)
415:       DOUBLE PRECISION, ALLOCATABLE :: MSBE(:)414:       DOUBLE PRECISION, ALLOCATABLE :: MSBE(:)
416:       DOUBLE PRECISION, ALLOCATABLE :: GAUSSKK(:,:)415:       DOUBLE PRECISION, ALLOCATABLE :: GAUSSKK(:,:)
417:       DOUBLE PRECISION, ALLOCATABLE :: GAUSSEE(:), GKSMALL(:)416:       DOUBLE PRECISION, ALLOCATABLE :: GAUSSEE(:), GKSMALL(:)
418:       DOUBLE PRECISION, ALLOCATABLE :: CSMIMAGES(:)417:       DOUBLE PRECISION, ALLOCATABLE :: CSMIMAGES(:)
419:       DOUBLE PRECISION, ALLOCATABLE :: QMINPCSMAV(:,:), CSMAV(:), QMINAV(:), PTGP(:,:,:), PTGPGUIDE(:,:,:)418:       DOUBLE PRECISION, ALLOCATABLE :: QMINPCSMAV(:,:), CSMAV(:), QMINAV(:), PTGP(:,:,:), PTGPGUIDE(:,:,:)
420:       DOUBLE PRECISION, ALLOCATABLE :: LOWESTC(:)419:       DOUBLE PRECISION, ALLOCATABLE :: LOWESTC(:)


r31119/finalio.f90 2016-09-06 15:30:38.294436231 +0100 r31118/finalio.f90 2016-09-06 15:30:41.506479102 +0100
315:                 I20,' function calls. Point group ',A4, &315:                 I20,' function calls. Point group ',A4, &
316:                 ' of order ', I3,'.')316:                 ' of order ', I3,'.')
317:         ELSE ! Print stuff for HSA in OPTIM's min.data format317:         ELSE ! Print stuff for HSA in OPTIM's min.data format
318:            !318:            !
319:            NUM_ZERO_EVS=6319:            NUM_ZERO_EVS=6
320:            IF (NATOMS.EQ.2) NUM_ZERO_EVS=5320:            IF (NATOMS.EQ.2) NUM_ZERO_EVS=5
321:            IF (MIEFT) NUM_ZERO_EVS=0321:            IF (MIEFT) NUM_ZERO_EVS=0
322:            !322:            !
323:            IF (ALLOCATED(HESS)) DEALLOCATE(HESS)323:            IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
324:            ALLOCATE(HESS(3*NATOMS,3*NATOMS))324:            ALLOCATE(HESS(3*NATOMS,3*NATOMS))
325:            IF(NSPECIES(0)>1) THEN ! ds656> this also sets ATMASS325:            IF(NSPECIES(0)>1) & ! ds656> this also sets ATMASS
326:               CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)326:                 CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)
327:               CALL SET_LABELS(QMINT(J1,1:NATOMS),MYNODE+1) 
328:            ENDIF 
329:            CALL POTENTIAL(QMINP(J1,1:3*NATOMS), EVALS, EDUMMY, &327:            CALL POTENTIAL(QMINP(J1,1:3*NATOMS), EVALS, EDUMMY, &
330:                 .TRUE., .TRUE.) ! EVALS is dummy for gradient here328:                 .TRUE., .TRUE.) ! EVALS is dummy for gradient here
331:            IF(ABS(EDUMMY-QMIN(J1)) > ECONV) THEN329:            IF(ABS(EDUMMY-QMIN(J1)) > ECONV) THEN
332:               WRITE(MYUNIT,'(A)') &330:               WRITE(MYUNIT,'(A)') &
333:                    'finalio> WARNING: Unexpected energy value! Configuration changed on final requench?'331:                    'finalio> WARNING: Unexpected energy value! Configuration changed on final requench?'
334:            ENDIF332:            ENDIF
335:            CALL MASSWT() ! Weight the Hessian by mass. 333:            CALL MASSWT() ! Weight the Hessian by mass. 
336:            !334:            !
337:            WRITE(MYUNIT, '(A)') &335:            WRITE(MYUNIT, '(A)') &
338:                 'finalio> Using DSYEV to calculate Hessian eigenvalues'336:                 'finalio> Using DSYEV to calculate Hessian eigenvalues'
732:         !730:         !
733:         ! this writes 'lowest' in xyz (Xmakemol) format731:         ! this writes 'lowest' in xyz (Xmakemol) format
734:         !732:         !
735:         WRITE(MYUNIT,'(A,I6,A)') ' in finalio BLN block MYUNIT2=',MYUNIT2733:         WRITE(MYUNIT,'(A,I6,A)') ' in finalio BLN block MYUNIT2=',MYUNIT2
736:         DO J2=1,NATOMS734:         DO J2=1,NATOMS
737:            WRITE(MYUNIT2,'(2A1,1X,3F20.10)') BEADLETTER(J2),'L',(QMINP(J1,3*(J2-1)+J3),J3=1,3)735:            WRITE(MYUNIT2,'(2A1,1X,3F20.10)') BEADLETTER(J2),'L',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
738:         ENDDO736:         ENDDO
739:      ELSE IF(QALCST.OR.MULTIPERMT.OR.MLJT.OR.MGUPTAT.OR.MSCT) THEN737:      ELSE IF(QALCST.OR.MULTIPERMT.OR.MLJT.OR.MGUPTAT.OR.MSCT) THEN
740:         !738:         !
741:         CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)739:         CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)
742:         CALL SET_LABELS(QMINT(J1,1:NATOMS),MYNODE+1) 
743:         IF(SPECMASST) THEN 740:         IF(SPECMASST) THEN 
744:            EDUMMY = SUM(ATMASS)741:            EDUMMY = SUM(ATMASS)
745:            WRITE(MYUNIT,'(A,I6,A,E20.10)') &742:            WRITE(MYUNIT,'(A,I6,A,E20.10)') &
746:                 'finalio> Min ',J1,' has total mass ', EDUMMY743:                 'finalio> Min ',J1,' has total mass ', EDUMMY
747:         ENDIF744:         ENDIF
748:         !745:         !
749:         IF(STRESST) CALL CALC_STRESS(QMINP(J1,1:3*NATOMS),EDUMMY)746:         IF(STRESST) CALL CALC_STRESS(QMINP(J1,1:3*NATOMS),EDUMMY)
750:         !747:         !
751:         DO J2=1,NSPECIES(0) !ds656> Should not exceed 10748:         DO J2=1,NSPECIES(0) !ds656> Should not exceed 10
752:            IF(SPECLABELST) THEN749:            IF(SPECLABELST) THEN


r31119/finalq.f 2016-09-06 15:30:38.546439595 +0100 r31118/finalq.f 2016-09-06 15:30:41.810483160 +0100
 70:       WRITE(MYUNIT,'(A)') 'NOTE: these may NOT match the other output files - see below for a sorted list of Lowest minima' 70:       WRITE(MYUNIT,'(A)') 'NOTE: these may NOT match the other output files - see below for a sorted list of Lowest minima'
 71:       DO J1=1,NSAVE 71:       DO J1=1,NSAVE
 72:          IF (QMIN(J1).LT.1.0D10) THEN 72:          IF (QMIN(J1).LT.1.0D10) THEN
 73:             NATOMS=QMINNATOMS(J1) 73:             NATOMS=QMINNATOMS(J1)
 74:             DO J2=1,3*NATOMS 74:             DO J2=1,3*NATOMS
 75:                COORDS(J2,NP)=QMINP(J1,J2) 75:                COORDS(J2,NP)=QMINP(J1,J2)
 76:             ENDDO 76:             ENDDO
 77:              77:             
 78:             !ds656> 78:             !ds656>
 79:             !write(MYUNIT,*) (QMINT(J1,J2),J2=1,NATOMS) 79:             !write(MYUNIT,*) (QMINT(J1,J2),J2=1,NATOMS)
 80:             IF(NSPECIES(0)>1) THEN 80:             IF(NSPECIES(0)>1) CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)
 81:                CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1) 
 82:                CALL SET_LABELS(QMINT(J1,1:NATOMS),NP) 
 83:             ENDIF 
 84:             !<ds656 81:             !<ds656
 85:              82:             
 86:             NQ(NP)=NQ(NP)+1 83:             NQ(NP)=NQ(NP)+1
 87:             IF(CHRMMT.AND.ACESOLV) NCHENCALLS=ACEUPSTEP-1 84:             IF(CHRMMT.AND.ACESOLV) NCHENCALLS=ACEUPSTEP-1
 88: !fh301>{{{ 85: !fh301>{{{
 89:             STRUC = NQ(1) 86:             STRUC = NQ(1)
 90:             COUT=.TRUE. 87:             COUT=.TRUE.
 91: !fh301>}}} 88: !fh301>}}}
 92: ! hk286 89: ! hk286
 93:  90: 


r31119/initialization.f90 2016-09-06 15:30:38.806443065 +0100 r31118/initialization.f90 2016-09-06 15:30:42.082486791 +0100
321:    NNLIST_MAX = 12321:    NNLIST_MAX = 12
322:    ALLOCATE(NNLISTS(1:NATOMS,0:NNLIST_MAX))322:    ALLOCATE(NNLISTS(1:NATOMS,0:NNLIST_MAX))
323:    IF(QALCS_SURFT) THEN323:    IF(QALCS_SURFT) THEN
324:       NVSITES_MAX = 6*NATOMS !6*NINT(DBLE(NATOMS)**(2.0D0/3.0D0))324:       NVSITES_MAX = 6*NATOMS !6*NINT(DBLE(NATOMS)**(2.0D0/3.0D0))
325:       ALLOCATE(VSITES(NVSITES_MAX,3))325:       ALLOCATE(VSITES(NVSITES_MAX,3))
326:    ENDIF326:    ENDIF
327: ENDIF327: ENDIF
328: !328: !
329: ! set all atoms as mobile and partition by type329: ! set all atoms as mobile and partition by type
330: CALL RESET_ATOMLISTS(1) 330: CALL RESET_ATOMLISTS(1) 
331: DO J1=1,NPAR 
332:    CALL RESET_LABELS(J1) 
333: ENDDO 
334: !331: !
335: ! ds656> If atomic labels are not to be permuted, we set QMINT 332: ! ds656> If atomic labels are not to be permuted, we set QMINT 
336: ! here and it will remain fixed throughout.333: ! here and it will remain fixed throughout.
337: IF(.NOT. QALCST) THEN334: IF(.NOT. QALCST) THEN
338:    DO J1=1,NSAVE335:    DO J1=1,NSAVE
339:       QMINT(J1,1:NATOMS) = INVATOMLISTS(1:NATOMS,1)336:       QMINT(J1,1:NATOMS) = INVATOMLISTS(1:NATOMS,1)
340:    END DO337:    END DO
341: ENDIF338: ENDIF
342: !339: !
343: ! ds656> If species labels have been specified, the number of 340: ! ds656> If species labels have been specified, the number of 


r31119/keywords.f 2016-09-06 15:30:39.070446588 +0100 r31118/keywords.f 2016-09-06 15:30:42.362490528 +0100
219:          DUMPINTEOS=.FALSE.219:          DUMPINTEOS=.FALSE.
220:          DUMPINTXYZFREQ=100220:          DUMPINTXYZFREQ=100
221:          DUMPINTEOSFREQ=100221:          DUMPINTEOSFREQ=100
222:          KINT=0.0D0222:          KINT=0.0D0
223:          QCIAMBERT=.FALSE.223:          QCIAMBERT=.FALSE.
224: 224: 
225: ! ns644>    Adding the TWIST keyword:225: ! ns644>    Adding the TWIST keyword:
226:       TWISTT=.FALSE.226:       TWISTT=.FALSE.
227:       NTWISTGROUPS=0227:       NTWISTGROUPS=0
228: 228: 
229: ! ds656>    COORDS and LABELS allocated here229: 
230:       ALLOCATE(FIXSTEP(1),FIXTEMP(1),FIXBOTH(1),TEMP(1),ACCRAT(1),STEP(1),ASTEP(1),OSTEP(1),BLOCK(1),NT(1),NQ(1),EPREV(1),230:       ALLOCATE(FIXSTEP(1),FIXTEMP(1),FIXBOTH(1),TEMP(1),ACCRAT(1),STEP(1),ASTEP(1),OSTEP(1),BLOCK(1),NT(1),NQ(1),EPREV(1),
231:      @         JUMPMOVE(1),JUMPINT(1),JDUMP(1),COORDS(3*NATOMSALLOC,1),LABELS(NATOMSALLOC,1), COORDSO(3*NATOMSALLOC,1), LABELSO(NATOMSALLOC,1), 231:      @         JUMPMOVE(1),JUMPINT(1),JDUMP(1),COORDS(3*NATOMSALLOC,1),COORDSO(3*NATOMSALLOC,1),VAT(NATOMSALLOC,1),
232:      @         VAT(NATOMSALLOC,1), VATO(NATOMSALLOC,1),  232:      @         VATO(NATOMSALLOC,1),  
233:      @         JUMPTO(1),SHELLMOVES(1),PTGROUP(1),NSURFMOVES(1),NCORE(1))233:      @         JUMPTO(1),SHELLMOVES(1),PTGROUP(1),NSURFMOVES(1),NCORE(1))
234:       DO JP=1,1234:       DO JP=1,1
235:          EPREV(JP)=1.0D100235:          EPREV(JP)=1.0D100
236:          FIXSTEP(JP)=.FALSE.236:          FIXSTEP(JP)=.FALSE.
237:          FIXTEMP(JP)=.FALSE.237:          FIXTEMP(JP)=.FALSE.
238:          FIXBOTH(JP)=.FALSE.238:          FIXBOTH(JP)=.FALSE.
239:          TEMP(JP)=0.3D0239:          TEMP(JP)=0.3D0
240:          ACCRAT(JP)=0.5D0240:          ACCRAT(JP)=0.5D0
241:          STEP(JP)=0.3D0241:          STEP(JP)=0.3D0
242:          ASTEP(JP)=0.3D0242:          ASTEP(JP)=0.3D0
4623: 4623: 
4624:          CALL create_mtargets_array(TARGETS)4624:          CALL create_mtargets_array(TARGETS)
4625: 4625: 
4626:          DEALLOCATE(TARGETS)4626:          DEALLOCATE(TARGETS)
4627: 4627: 
4628: !4628: !
4629: !  MPI keyword4629: !  MPI keyword
4630: !4630: !
4631:       ELSE IF (WORD.EQ.'MPI') THEN4631:       ELSE IF (WORD.EQ.'MPI') THEN
4632:          MPIT=.TRUE.4632:          MPIT=.TRUE.
4633:          DEALLOCATE(FIXSTEP,FIXTEMP,FIXBOTH,TEMP,ACCRAT,STEP,ASTEP,OSTEP,BLOCK,NT,JUMPMOVE,JUMPINT,JDUMP,COORDS,LABELS,NQ,4633:          DEALLOCATE(FIXSTEP,FIXTEMP,FIXBOTH,TEMP,ACCRAT,STEP,ASTEP,OSTEP,BLOCK,NT,JUMPMOVE,JUMPINT,JDUMP,COORDS,NQ,
4634:      @              JUMPTO,EPREV,COORDSO,LABELSO,VAT,VATO,SHELLMOVES,PTGROUP,NSURFMOVES,NCORE)4634:      @              JUMPTO,EPREV,COORDSO,VAT,VATO,SHELLMOVES,PTGROUP,NSURFMOVES,NCORE)
4635:          ALLOCATE(FIXSTEP(NPAR),FIXTEMP(NPAR),FIXBOTH(NPAR),TEMP(NPAR),ACCRAT(NPAR),STEP(NPAR),ASTEP(NPAR),OSTEP(NPAR),4635:          ALLOCATE(FIXSTEP(NPAR),FIXTEMP(NPAR),FIXBOTH(NPAR),TEMP(NPAR),ACCRAT(NPAR),STEP(NPAR),ASTEP(NPAR),OSTEP(NPAR),
4636:      @         BLOCK(NPAR),NT(NPAR),JUMPMOVE(NPAR),JUMPINT(NPAR),JDUMP(NPAR),COORDS(3*NATOMSALLOC,NPAR),LABELS(NATOMSALLOC,NPAR),NQ(NPAR),4636:      @         BLOCK(NPAR),NT(NPAR),JUMPMOVE(NPAR),JUMPINT(NPAR),JDUMP(NPAR),COORDS(3*NATOMSALLOC,NPAR),NQ(NPAR),
4637:      @         JUMPTO(NPAR),EPREV(NPAR),4637:      @         JUMPTO(NPAR),EPREV(NPAR),
4638:      @         COORDSO(3*NATOMSALLOC,NPAR),LABELSO(NATOMSALLOC,NPAR), VAT(NATOMSALLOC,NPAR),VATO(NATOMSALLOC,NPAR))4638:      @         COORDSO(3*NATOMSALLOC,NPAR),VAT(NATOMSALLOC,NPAR),VATO(NATOMSALLOC,NPAR))
4639:          ALLOCATE(SHELLMOVES(NPAR))4639:          ALLOCATE(SHELLMOVES(NPAR))
4640:          ALLOCATE(PTGROUP(NPAR))4640:          ALLOCATE(PTGROUP(NPAR))
4641:          ALLOCATE(NSURFMOVES(NPAR))4641:          ALLOCATE(NSURFMOVES(NPAR))
4642:          ALLOCATE(NCORE(NPAR))4642:          ALLOCATE(NCORE(NPAR))
4643:          DO JP=1,NPAR4643:          DO JP=1,NPAR
4644:             EPREV(JP)=1.0D1004644:             EPREV(JP)=1.0D100
4645:             FIXSTEP(JP)=.FALSE.4645:             FIXSTEP(JP)=.FALSE.
4646:             FIXTEMP(JP)=.FALSE.4646:             FIXTEMP(JP)=.FALSE.
4647:             FIXBOTH(JP)=.FALSE.4647:             FIXBOTH(JP)=.FALSE.
4648:             TEMP(JP)=0.3D04648:             TEMP(JP)=0.3D0
4929:             WRITE(MYUNIT,*) PAIRDIST(J1,:)4929:             WRITE(MYUNIT,*) PAIRDIST(J1,:)
4930:          ENDDO4930:          ENDDO
4931: 4931: 
4932: !4932: !
4933: !  PARALLEL must come before STEP and ACCRAT4933: !  PARALLEL must come before STEP and ACCRAT
4934: !  This keyword is for the serial parallel implementation - now obsolete.4934: !  This keyword is for the serial parallel implementation - now obsolete.
4935: !4935: !
4936:       ELSE IF (WORD.EQ.'PARALLEL') THEN4936:       ELSE IF (WORD.EQ.'PARALLEL') THEN
4937:          PARALLELT=.TRUE.4937:          PARALLELT=.TRUE.
4938:          CALL READI(NPAR)4938:          CALL READI(NPAR)
4939:          DEALLOCATE(FIXSTEP,FIXTEMP,FIXBOTH,TEMP,ACCRAT,STEP,ASTEP,OSTEP,BLOCK,NT,JUMPMOVE,JUMPINT,JDUMP,COORDS,LABELS,NQ,4939:          DEALLOCATE(FIXSTEP,FIXTEMP,FIXBOTH,TEMP,ACCRAT,STEP,ASTEP,OSTEP,BLOCK,NT,JUMPMOVE,JUMPINT,JDUMP,COORDS,NQ,
4940:      @              JUMPTO,EPREV,COORDSO,LABELSO,VAT,VATO,SHELLMOVES,PTGROUP,NSURFMOVES,NCORE) 4940:      @              JUMPTO,EPREV,COORDSO,VAT,VATO,SHELLMOVES,PTGROUP,NSURFMOVES,NCORE) 
4941:          ALLOCATE(FIXSTEP(NPAR),FIXTEMP(NPAR),FIXBOTH(NPAR),TEMP(NPAR),ACCRAT(NPAR),STEP(NPAR),ASTEP(NPAR),OSTEP(NPAR), 4941:          ALLOCATE(FIXSTEP(NPAR),FIXTEMP(NPAR),FIXBOTH(NPAR),TEMP(NPAR),ACCRAT(NPAR),STEP(NPAR),ASTEP(NPAR),OSTEP(NPAR), 
4942:      @         BLOCK(NPAR),NT(NPAR),JUMPMOVE(NPAR),JUMPINT(NPAR),JDUMP(NPAR),NQ(NPAR),JUMPTO(NPAR),4942:      @         BLOCK(NPAR),NT(NPAR),JUMPMOVE(NPAR),JUMPINT(NPAR),JDUMP(NPAR),NQ(NPAR),JUMPTO(NPAR),
4943:      &         COORDS(3*NATOMSALLOC,NPAR), LABELS(NATOMSALLOC,NPAR),4943:      &         COORDS(3*NATOMSALLOC,NPAR),
4944:      @         COORDSO(3*NATOMSALLOC,NPAR), LABELSO(NATOMSALLOC,NPAR), VAT(NATOMSALLOC,NPAR),VATO(NATOMSALLOC,NPAR),EPREV(NPAR),SHELLMOVES(NPAR),PTGROUP(NPAR),4944:      @         COORDSO(3*NATOMSALLOC,NPAR),VAT(NATOMSALLOC,NPAR),VATO(NATOMSALLOC,NPAR),EPREV(NPAR),SHELLMOVES(NPAR),PTGROUP(NPAR),
4945:      @         NSURFMOVES(NPAR),NCORE(NPAR),REPLOW(NPAR))4945:      @         NSURFMOVES(NPAR),NCORE(NPAR),REPLOW(NPAR))
4946:          NATOMS=NATOMS/NPAR4946:          NATOMS=NATOMS/NPAR
4947:          DO JP=1,NPAR4947:          DO JP=1,NPAR
4948:             FIXSTEP(JP)=.FALSE.4948:             FIXSTEP(JP)=.FALSE.
4949:             FIXTEMP(JP)=.FALSE.4949:             FIXTEMP(JP)=.FALSE.
4950:             FIXBOTH(JP)=.FALSE.4950:             FIXBOTH(JP)=.FALSE.
4951:             TEMP(JP)=0.3D04951:             TEMP(JP)=0.3D0
4952:             ACCRAT(JP)=0.5D04952:             ACCRAT(JP)=0.5D0
4953:             STEP(JP)=0.3D04953:             STEP(JP)=0.3D0
4954:             ASTEP(JP)=0.3D04954:             ASTEP(JP)=0.3D0


r31119/mc.F 2016-09-06 15:30:39.350450326 +0100 r31118/mc.F 2016-09-06 15:30:42.770495977 +0100
343:       IF (DEBUG .AND. PERCOLATET) PERCCOMPMARKOV = COMPON343:       IF (DEBUG .AND. PERCOLATET) PERCCOMPMARKOV = COMPON
344: ! jwrm2> end344: ! jwrm2> end
345:       EPPREV(JP)=0.0D0345:       EPPREV(JP)=0.0D0
346:       ELASTSYM(JP)=0.0D0346:       ELASTSYM(JP)=0.0D0
347:       IF (.NOT.RESTORET) EBEST(JP)=POTEL347:       IF (.NOT.RESTORET) EBEST(JP)=POTEL
348:       BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)348:       BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
349:       JBEST(JP)=0349:       JBEST(JP)=0
350:       RMIN=POTEL350:       RMIN=POTEL
351:       RCOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,1)351:       RCOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,1)
352:       COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)352:       COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
353:       LABELSO(1:NATOMS,JP)=LABELS(1:NATOMS,JP) ! <ds656 
354:       VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)353:       VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)
355:       EPSSPHERE=EPSSAVE354:       EPSSPHERE=EPSSAVE
356: 355: 
357: ! Initialisation 356: ! Initialisation 
358: 357: 
359:       IF (PTTMIN < 1.0D-6 ) PTTMIN = 1.0D-6 ! to avoid devision by zero358:       IF (PTTMIN < 1.0D-6 ) PTTMIN = 1.0D-6 ! to avoid devision by zero
360:       CTE=(LOG(PTTMAX/PTTMIN))/(NPAR-1)359:       CTE=(LOG(PTTMAX/PTTMIN))/(NPAR-1)
361:       CTE=EXP(CTE)360:       CTE=EXP(CTE)
362: 361: 
363:       DO I=0, NPAR-1362:       DO I=0, NPAR-1
466:          IF (DEBUG .AND. PERCOLATET) PERCCOMPMARKOV = COMPON465:          IF (DEBUG .AND. PERCOLATET) PERCCOMPMARKOV = COMPON
467: ! jwrm2> end466: ! jwrm2> end
468:          EPPREV(JP)=0.0D0467:          EPPREV(JP)=0.0D0
469:          ELASTSYM(JP)=0.0D0468:          ELASTSYM(JP)=0.0D0
470:          IF (.NOT.RESTORET) EBEST(JP)=POTEL469:          IF (.NOT.RESTORET) EBEST(JP)=POTEL
471:          BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)470:          BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
472:          JBEST(JP)=0471:          JBEST(JP)=0
473:          RMIN=POTEL472:          RMIN=POTEL
474:          RCOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,1)473:          RCOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,1)
475:          COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)474:          COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
476:          LABELSO(1:NATOMS,JP)=LABELS(1:NATOMS,JP) ! <ds656 
477:          VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)475:          VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)
478: !hk286 - otherwise NEWRESTART is always called when the keyword RESTORE is also used476: !hk286 - otherwise NEWRESTART is always called when the keyword RESTORE is also used
479:          IF (.NOT.RESTORET) THEN477:          IF (.NOT.RESTORET) THEN
480:             EBEST(JP)=POTEL478:             EBEST(JP)=POTEL
481:             JBEST(JP)=0479:             JBEST(JP)=0
482:             BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)480:             BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
483:          ENDIF481:          ENDIF
484: !482: !
485: ! Set up data for relxation blocks with grand canonical BH483: ! Set up data for relxation blocks with grand canonical BH
486: !484: !
1267:                   IF (CHANGED) THEN1265:                   IF (CHANGED) THEN
1268:                      CALL QUENCH(.FALSE.,JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)1266:                      CALL QUENCH(.FALSE.,JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
1269:                      IF (RMS.GT.BQMAX) THEN1267:                      IF (RMS.GT.BQMAX) THEN
1270:                         WRITE(MYUNIT,'(A)') 'mc> Quench from new size geometry failed - reject GC step attempt'1268:                         WRITE(MYUNIT,'(A)') 'mc> Quench from new size geometry failed - reject GC step attempt'
1271:                         WRITE(MYUNIT,'(A,2G20.10)') 'mc> RMS and BQMAX=',RMS,BQMAX1269:                         WRITE(MYUNIT,'(A,2G20.10)') 'mc> RMS and BQMAX=',RMS,BQMAX
1272:                         CHANGED=.FALSE.1270:                         CHANGED=.FALSE.
1273:                         IF (NATOMS.GT.GCNATOMSBEST(JP)) SYMFCTR=SYMFCTR-LOG(1.0D0*NATOMS)1271:                         IF (NATOMS.GT.GCNATOMSBEST(JP)) SYMFCTR=SYMFCTR-LOG(1.0D0*NATOMS)
1274:                         IF (NATOMS.LT.GCNATOMSBEST(JP)) SYMFCTR=SYMFCTR+LOG(1.0D0*GCNATOMSBEST(JP))1272:                         IF (NATOMS.LT.GCNATOMSBEST(JP)) SYMFCTR=SYMFCTR+LOG(1.0D0*GCNATOMSBEST(JP))
1275:                         NATOMS=GCNATOMSBEST(JP)1273:                         NATOMS=GCNATOMSBEST(JP)
1276:                         COORDS(1:3*NATOMS,JP)=COORDSO(1:3*NATOMS,JP)1274:                         COORDS(1:3*NATOMS,JP)=COORDSO(1:3*NATOMS,JP)
1277:                         LABELS(1:NATOMS,JP)=LABELSO(1:NATOMS,JP) ! <ds656 
1278:                         VAT(1:NATOMS,JP)=VATO(1:NATOMS,JP)1275:                         VAT(1:NATOMS,JP)=VATO(1:NATOMS,JP)
1279:                         POTEL=EPREV(JP)1276:                         POTEL=EPREV(JP)
1280:                         WRITE(MYUNIT,'(A,I6)') 'mc> Number of atoms reset to ',NATOMS1277:                         WRITE(MYUNIT,'(A,I6)') 'mc> Number of atoms reset to ',NATOMS
1281:                      ENDIF1278:                      ENDIF
1282: 1279: 
1283:                      IF (PERCOLATET) THEN1280:                      IF (PERCOLATET) THEN
1284:                         IF (.NOT. PERCT) THEN1281:                         IF (.NOT. PERCT) THEN
1285:                            WRITE(MYUNIT,'(A)') 'mc> Disconnected atoms, rejecting step. '1282:                            WRITE(MYUNIT,'(A)') 'mc> Disconnected atoms, rejecting step. '
1286:                            CHANGED=.FALSE.1283:                            CHANGED=.FALSE.
1287:                            IF (NATOMS.GT.GCNATOMSBEST(JP)) SYMFCTR=SYMFCTR-LOG(1.0D0*NATOMS)1284:                            IF (NATOMS.GT.GCNATOMSBEST(JP)) SYMFCTR=SYMFCTR-LOG(1.0D0*NATOMS)
1288:                            IF (NATOMS.LT.GCNATOMSBEST(JP)) SYMFCTR=SYMFCTR+LOG(1.0D0*GCNATOMSBEST(JP))1285:                            IF (NATOMS.LT.GCNATOMSBEST(JP)) SYMFCTR=SYMFCTR+LOG(1.0D0*GCNATOMSBEST(JP))
1289:                            NATOMS=GCNATOMSBEST(JP)1286:                            NATOMS=GCNATOMSBEST(JP)
1290:                            COORDS(1:3*NATOMS,JP)=COORDSO(1:3*NATOMS,JP)1287:                            COORDS(1:3*NATOMS,JP)=COORDSO(1:3*NATOMS,JP)
1291:                            LABELS(1:NATOMS,JP)=LABELSO(1:NATOMS,JP) ! <ds656 
1292:                            VAT(1:NATOMS,JP)=VATO(1:NATOMS,JP)1288:                            VAT(1:NATOMS,JP)=VATO(1:NATOMS,JP)
1293:                            POTEL=EPREV(JP)1289:                            POTEL=EPREV(JP)
1294:                            WRITE(MYUNIT,'(A,I6)') 'mc> Number of atoms reset to ',NATOMS1290:                            WRITE(MYUNIT,'(A,I6)') 'mc> Number of atoms reset to ',NATOMS
1295:                         ENDIF1291:                         ENDIF
1296:                      ENDIF1292:                      ENDIF
1297:                      IF (MOD(J1-1,PRTFRQ).EQ.0) WRITE(MYUNIT,'(A,G20.10)') 'mc> New quench energy=',POTEL1293:                      IF (MOD(J1-1,PRTFRQ).EQ.0) WRITE(MYUNIT,'(A,G20.10)') 'mc> New quench energy=',POTEL
1298:                   ENDIF1294:                   ENDIF
1299: 1295: 
1300:                   IF (CHANGED) THEN1296:                   IF (CHANGED) THEN
1301: !1297: !
1317:                      NSUCCESS(JP)=01313:                      NSUCCESS(JP)=0
1318:                      NFAIL(JP)=01314:                      NFAIL(JP)=0
1319:                      EBEST(JP)=POTEL ! this is communicated via common block MYPOT, which is in quench.f1315:                      EBEST(JP)=POTEL ! this is communicated via common block MYPOT, which is in quench.f
1320:                      BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)1316:                      BESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
1321:                      JBEST(JP)=J11317:                      JBEST(JP)=J1
1322:                      EPREV(JP)=POTEL1318:                      EPREV(JP)=POTEL
1323:                      EPPREV(JP)=0.0D01319:                      EPPREV(JP)=0.0D0
1324:                      NSYMREM=01320:                      NSYMREM=0
1325:                      IF (CENT.AND.(.NOT.SEEDT)) CALL CENTRE2(COORDS(1:3*NATOMS,JP))1321:                      IF (CENT.AND.(.NOT.SEEDT)) CALL CENTRE2(COORDS(1:3*NATOMS,JP))
1326:                      COORDSO(1:3*(NATOMS-NSEED),JP)=COORDS(1:3*(NATOMS-NSEED),JP)1322:                      COORDSO(1:3*(NATOMS-NSEED),JP)=COORDS(1:3*(NATOMS-NSEED),JP)
1327:                      LABELSO(1:(NATOMS-NSEED),JP)=LABELS(1:(NATOMS-NSEED),JP) ! <ds656 
1328:                      VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)1323:                      VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)
1329: 1324: 
1330:                      GCEBEST(JP)=POTEL1325:                      GCEBEST(JP)=POTEL
1331:                      GCJBEST(JP)=J11326:                      GCJBEST(JP)=J1
1332:                      GCBESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)1327:                      GCBESTCOORDS(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
1333:                      GCBESTVAT(1:NATOMS,JP)=VAT(1:NATOMS,JP)1328:                      GCBESTVAT(1:NATOMS,JP)=VAT(1:NATOMS,JP)
1334:                      GCNATOMSBEST(JP)=NATOMS1329:                      GCNATOMSBEST(JP)=NATOMS
1335:                   ENDIF1330:                   ENDIF
1336: !1331: !
1337: ! ALL OTHER STEP TAKING1332: ! ALL OTHER STEP TAKING
2070:                            HBONDGROUPS(HBONDMATCHN,:,:)=HBONDMAT2065:                            HBONDGROUPS(HBONDMATCHN,:,:)=HBONDMAT
2071: ! csw34> Dump an updates group*.mat file for the group2066: ! csw34> Dump an updates group*.mat file for the group
2072:                            WRITE(NHBONDGROUPSCHAR,*) HBONDMATCHN2067:                            WRITE(NHBONDGROUPSCHAR,*) HBONDMATCHN
2073:                            HBONDGROUPNAME='group'//TRIM(ADJUSTL(NHBONDGROUPSCHAR))2068:                            HBONDGROUPNAME='group'//TRIM(ADJUSTL(NHBONDGROUPSCHAR))
2074:                            CALL WRITEHBONDMATRIX(HBONDMAT(:,:),HBONDNRES,HBONDGROUPNAME)2069:                            CALL WRITEHBONDMATRIX(HBONDMAT(:,:),HBONDNRES,HBONDGROUPNAME)
2075:                         ENDIF2070:                         ENDIF
2076:                      ENDIF2071:                      ENDIF
2077:                   ENDIF2072:                   ENDIF
2078: ! END OF HBONDMATRIX2073: ! END OF HBONDMATRIX
2079:                   COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)2074:                   COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
2080:                   LABELSO(1:NATOMS,JP)=LABELS(1:NATOMS,JP) ! <ds656 
2081:                   VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)2075:                   VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)
2082:                ELSE2076:                ELSE
2083:                   NFAIL(JP)=NFAIL(JP)+12077:                   NFAIL(JP)=NFAIL(JP)+1
2084:                   REJSTREAK=REJSTREAK+12078:                   REJSTREAK=REJSTREAK+1
2085:                   MARKOVSTEP=.FALSE.2079:                   MARKOVSTEP=.FALSE.
2086:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)2080:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)
2087:                   IF (DEBUG) THEN2081:                   IF (DEBUG) THEN
2088:                      WRITE(MYUNIT,36) JP,RANDOM,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)2082:                      WRITE(MYUNIT,36) JP,RANDOM,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)
2089: 36                   FORMAT('JP,RAN,POTEL,EPREV,NSUC,NFAIL=',I2,3F15.7,2I6,' REJ')2083: 36                   FORMAT('JP,RAN,POTEL,EPREV,NSUC,NFAIL=',I2,3F15.7,2I6,' REJ')
2090:                   ENDIF2084:                   ENDIF
2571:             READ(UNT,*) (COORDS(J2,JP),J2=1,3*NATOMS)2565:             READ(UNT,*) (COORDS(J2,JP),J2=1,3*NATOMS)
2572: !2566: !
2573: !  Coordinates should be converged already, but we need to reset VAT and VATO.2567: !  Coordinates should be converged already, but we need to reset VAT and VATO.
2574: !2568: !
2575: !  next line should be uncommented if routine is made availabe to use with CHARMM2569: !  next line should be uncommented if routine is made availabe to use with CHARMM
2576: !            IF(CHRMMT.AND.ACESOLV) NCHENCALLS=ACEUPSTEP-12570: !            IF(CHRMMT.AND.ACESOLV) NCHENCALLS=ACEUPSTEP-1
2577:             CALL QUENCH(.FALSE.,JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)2571:             CALL QUENCH(.FALSE.,JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
2578: 2572: 
2579:             DO J2=1,NATOMS2573:             DO J2=1,NATOMS
2580:                VATO(J2,JP)=VAT(J2,JP)2574:                VATO(J2,JP)=VAT(J2,JP)
2581:                LABELSO(J2,JP)=LABELS(J2,JP) !< ds656 
2582:             ENDDO2575:             ENDDO
2583:             DO J2=1,3*NATOMS2576:             DO J2=1,3*NATOMS
2584:                COORDSO(J2,JP)=COORDS(J2,JP)2577:                COORDSO(J2,JP)=COORDS(J2,JP)
2585:             ENDDO2578:             ENDDO
2586:             IF (DEBUG) THEN2579:             IF (DEBUG) THEN
2587:                WRITE(MYUNIT,'(A)') 'Jump coordinates:'2580:                WRITE(MYUNIT,'(A)') 'Jump coordinates:'
2588:                WRITE(MYUNIT,'(3F20.10)') (COORDS(J2,JP),J2=1,3*NATOMS)2581:                WRITE(MYUNIT,'(3F20.10)') (COORDS(J2,JP),J2=1,3*NATOMS)
2589:             ENDIF2582:             ENDIF
2590:             DO J2=1,NATOMS*(NQ(JUMPTO(JP))-1)-NDUM*NATOMS2583:             DO J2=1,NATOMS*(NQ(JUMPTO(JP))-1)-NDUM*NATOMS
2591:                READ(UNT,*) DUMMY2584:                READ(UNT,*) DUMMY


r31119/MGupta.f90 2016-09-06 15:30:36.670414551 +0100 r31118/MGupta.f90 2016-09-06 15:30:39.890457533 +0100
 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,GRAD,POT,GRADT,HESST,STRESST)
 26:   ! 26:   !
 27:   USE COMMONS, ONLY : NATOMS, VT, NSPECIES, ATOMLISTS,& 27:   USE COMMONS, ONLY : NATOMS, VT, NSPECIES, ATOMLISTS,&
 28:        INVATOMLISTS, STRESS, MYNODE, LABELS 28:        INVATOMLISTS, STRESS
 29:   USE MODHESS, ONLY : HESS 29:   USE MODHESS, ONLY : HESS
 30:   USE POT_PARAMS, ONLY : MGUPTA_A, MGUPTA_XI, MGUPTA_P, & 30:   USE POT_PARAMS, ONLY : MGUPTA_A, MGUPTA_XI, MGUPTA_P, &
 31:        MGUPTA_Q, MGUPTA_R0, MGUPTA_M2Q, MGUPTA_XI_SQ, & 31:        MGUPTA_Q, MGUPTA_R0, MGUPTA_M2Q, MGUPTA_XI_SQ, &
 32:        MGUPTA_MP_DIVBY_R0, MGUPTA_M2Q_DIVBY_R0 32:        MGUPTA_MP_DIVBY_R0, MGUPTA_M2Q_DIVBY_R0
 33:   ! 33:   !
 34:   IMPLICIT NONE 34:   IMPLICIT NONE
 35:   ! 35:   !
 36:   ! Passed variables ----------------------------------------- 36:   ! Passed variables -----------------------------------------
 37:   LOGICAL, INTENT(IN) :: GRADT, HESST, STRESST 37:   LOGICAL, INTENT(IN) :: GRADT, HESST, STRESST
 38:   DOUBLE PRECISION, INTENT(IN) :: X(3*NATOMS) 38:   DOUBLE PRECISION, INTENT(IN) :: X(3*NATOMS)
 39:   DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), POT 39:   DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), POT
 40:   ! 40:   !
 41:   ! Local variables ----------------------------------------- 41:   ! Local variables -----------------------------------------
 42:   INTEGER :: I,J,J1,J2,J3,J13,J23,J33,T1,T2,T3,IP 42:   LOGICAL :: SELF
  43:   INTEGER :: I,J,J1,J2,J3,J13,J23,J33,T1,T2,T3,G1,G2,G3,&
  44:        K1,K2,K3,LI1,UI1,LI2,UI2,G2START
 43:   DOUBLE PRECISION :: DIST, DX(3), RHO(NATOMS), DUMMY, & 45:   DOUBLE PRECISION :: DIST, DX(3), RHO(NATOMS), DUMMY, &
 44:        RIJ(NATOMS,NATOMS), DUDR, D2UDR2, IGRAD(3), & 46:        RIJ(NATOMS,NATOMS), DUDR, D2UDR2, IGRAD(3), &
 45:        REP_IJ(NATOMS,NATOMS), ATT_IJ(NATOMS,NATOMS), & 47:        REP_IJ(NATOMS,NATOMS), ATT_IJ(NATOMS,NATOMS), &
 46:        REP_DUM1, REP_DUM2, REP_DUM3, ATT_DUM1, ATT_DUM2, & 48:        REP_DUM1, REP_DUM2, REP_DUM3, ATT_DUM1, ATT_DUM2, &
 47:        ATT_DUM3, RHO_DUM, RHO4HESS(NATOMS) 49:        ATT_DUM3, RHO_DUM, RHO4HESS(NATOMS)
 48:   ! 50:   !
 49:   ! Do some initialisations 51:   ! Do some initialisations
 50:   RIJ(1:NATOMS,1:NATOMS)=0.0D0 52:   RIJ(1:NATOMS,1:NATOMS)=0.0D0
 51:   REP_IJ(1:NATOMS,1:NATOMS)=0.0D0 53:   REP_IJ(1:NATOMS,1:NATOMS)=0.0D0
 52:   ATT_IJ(1:NATOMS,1:NATOMS)=0.0D0 54:   ATT_IJ(1:NATOMS,1:NATOMS)=0.0D0
 53:   RHO(1:NATOMS)=0.0D0 55:   RHO(1:NATOMS)=0.0D0
 54:   RHO4HESS(1:NATOMS)=0.0D0 56:   RHO4HESS(1:NATOMS)=0.0D0
 55:   VT(1:NATOMS) = 0.0D0 57:   VT(1:NATOMS) = 0.0D0
 56:   ! 58:   !
 57:   IP=MYNODE+1 ! Atoms labels from global LABELS(1:NATOMS,IP) 59:   ! Double loop over atom types
 58:   ! 60:   LT1: DO T1=1,NSPECIES(0)
 59:   DO J1=1,NATOMS-1 61:      LT2: DO T2=1,T1 ! T1,NSPECIES(0)
 60:      J13 = 3*(J1-1) 
 61:      T1=LABELS(J1,IP) 
 62:      DO J2=J1+1,NATOMS 
 63:         J23 = 3*(J2-1) 
 64:         T2=LABELS(J2,IP) 
 65:         ! 
 66:         DIST=0.0D0 
 67:         DO I=1,3 
 68:            DIST = DIST + (X(J23+I)-X(J13+I))**2 
 69:         ENDDO 
 70:         DIST = DSQRT(DIST) 
 71:         ! 
 72:         RIJ(J2,J1)=DIST ! store distance 
 73:         RIJ(J1,J2)=DIST ! impose symmetry 
 74:         ! 
 75:         DUMMY = MGUPTA_A(T1,T2)*& 
 76:              DEXP(MGUPTA_MP_DIVBY_R0(T1,T2)*DIST + & 
 77:              MGUPTA_P(T1,T2)) 
 78:         REP_IJ(J1,J2) = DUMMY 
 79:         REP_IJ(J2,J1) = DUMMY 
 80:         ! 
 81:         VT(J1) = VT(J1) + DUMMY 
 82:         VT(J2) = VT(J2) + DUMMY 
 83:         ! 62:         !
 84:         DUMMY = MGUPTA_XI_SQ(T1,T2)*& 63:         ! Double loop over mobile and frozen atom groups
 85:              DEXP(MGUPTA_M2Q_DIVBY_R0(T1,T2)*DIST - & 64:         LG1: DO G1 = 1,2 ! 1 -> mobile; 2 -> frozen
 86:              MGUPTA_M2Q(T1,T2)) 
 87:         ATT_IJ(J1,J2) = DUMMY 
 88:         ATT_IJ(J2,J1) = DUMMY 
 89:         ! 
 90:         ! Accumulate many-body potential term: 
 91:         RHO(J1) = RHO(J1) + DUMMY 
 92:         RHO(J2) = RHO(J2) + DUMMY 
 93:         ! 
 94:         IF (GRADT.OR.HESST.OR.STRESST) THEN 
 95:            ! 65:            !
 96:            ! Precompute derivatives of rep and att here 66:            IF(T1 == T2) THEN
 97:            REP_IJ(J1,J2) = MGUPTA_MP_DIVBY_R0(T1,T2)*& 67:               G2START = G1
 98:                 REP_IJ(J1,J2)*2.0D0  ! need the 2 here! 68:            ELSE
 99:            REP_IJ(J2,J1) = REP_IJ(J1,J2) 69:               G2START = 1
100:            ATT_IJ(J1,J2) = MGUPTA_M2Q_DIVBY_R0(T1,T2)*& 70:            ENDIF
101:                 ATT_IJ(J1,J2) 
102:            ATT_IJ(J2,J1) = ATT_IJ(J1,J2) 
103:            ! 71:            !
104:         ENDIF 72:            LG2: DO G2 = G2START,2 ! include frozen-frozen interactions
  73:               !
  74:               LI1=1
  75:               UI1=ATOMLISTS(T1,G1,0)
  76:               !
  77:               IF(T1==T2.AND.G1==G2) THEN
  78:                  SELF = .TRUE.
  79:                  UI1 = UI1 - 1 ! Exclude self-interaction
  80:               ELSE
  81:                  SELF = .FALSE.
  82:               ENDIF
  83:               !
  84:               LA1: DO K1=LI1,UI1
  85:                  !
  86:                  J1 = ATOMLISTS(T1,G1,K1) ! actual atom 1 index
  87:                  J13 = 3*(J1-1)
  88:                  !
  89:                  ! Atom loop bounds depend on type and group.
  90:                  ! Can this IF block be moved outside of loop LA1?
  91:                  IF(SELF) THEN
  92:                     LI2=K1+1
  93:                     UI2=UI1+1
  94:                  ELSE
  95:                     LI2=1
  96:                     UI2=ATOMLISTS(T2,G2,0)
  97:                  ENDIF
  98:                  !
  99:                  LA2: DO K2=LI2,UI2
 100:                     !
 101:                     J2 = ATOMLISTS(T2,G2,K2) ! actual atom 2 index
 102:                     J23 = 3*(J2-1)
 103:                     !
 104:                     !TESTCOUNT=TESTCOUNT+1
 105:                     !WRITE(*,*) "BGupta> TEST:", J1,J2,T1,T2
 106:                     !
 107:                     DIST=0.0D0
 108:                     DO I=1,3
 109:                        DIST = DIST + (X(J23+I)-X(J13+I))**2
 110:                     ENDDO
 111:                     DIST = DSQRT(DIST)
 112:                     !
 113:                     RIJ(J2,J1)=DIST ! store distance
 114:                     RIJ(J1,J2)=DIST ! impose symmetry
 115:                     !
 116:                     DUMMY = MGUPTA_A(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:                     !
 122:                     VT(J1) = VT(J1) + DUMMY
 123:                     VT(J2) = VT(J2) + DUMMY
 124:                     !
 125:                     DUMMY = MGUPTA_XI_SQ(T1,T2)*&
 126:                          DEXP(MGUPTA_M2Q_DIVBY_R0(T1,T2)*DIST - &
 127:                          MGUPTA_M2Q(T1,T2))
 128:                     ATT_IJ(J1,J2) = DUMMY
 129:                     ATT_IJ(J2,J1) = 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:                     !
 147:                  ENDDO LA2 ! Double loop over atoms
 148:               ENDDO LA1
 149:               !
 150:            ENDDO LG2 ! Span dynamic (1) and frozen (2) groups only
 151:         ENDDO LG1
105:         !152:         !
106:      ENDDO153:      ENDDO LT2 ! Double loop over atom types
107:   ENDDO154:   ENDDO LT1
108:   !155:   !
109:   ! Now store per-atom energies in array VT and sum over all atoms156:   ! Now store per-atom energies in array VT and sum over all atoms
110:   !157:   !
111:   POT=0.0D0158:   POT=0.0D0
112:   DO J1=1,NATOMS159:   DO G1=1,2 ! Groups 1 and 2
113:      RHO(J1)=DSQRT(RHO(J1)) ! square root density160:      DO T1=1,NSPECIES(0) ! Atom types
114:      VT(J1) = VT(J1) - RHO(J1) ! accummulate per-atom energy161:         DO J2=1,ATOMLISTS(T1,G1,0) 
115:      POT = POT + VT(J1)       ! accumulate potential energy162:            J1=ATOMLISTS(T1,G1,J2)
116:      IF(GRADT.OR.GRADT.OR.STRESST) THEN163:            RHO(J1)=DSQRT(RHO(J1)) ! square root density
117:         RHO(J1)=1.0D0/RHO(J1) 164:            VT(J1) = VT(J1) - RHO(J1) ! accummulate per-atom energy
118:         IF(HESST) RHO4HESS(J1) = 0.25D0*RHO(J1)**3165:            POT = POT + VT(J1)       ! accumulate potential energy
119:         RHO(J1)=0.5D0*RHO(J1)166:            IF(GRADT.OR.GRADT.OR.STRESST) THEN
120:      ENDIF167:               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:         ENDDO
 172:      ENDDO
121:   ENDDO173:   ENDDO
122:   !174:   !
123:   ! Calculate gradient terms, if required175:   ! Calculate gradient terms, if required
124:   !176:   !
125:   IF (GRADT.OR.HESST.OR.STRESST) THEN177:   IF (GRADT.OR.HESST.OR.STRESST) THEN
126:      ! 178:      ! 
127:      GRAD(:) = 0.0D0179:      GRAD(:) = 0.0D0
128:      IF(HESST) HESS(:,:) = 0.0D0180:      IF(HESST) HESS(:,:) = 0.0D0
129:      IF(STRESST) STRESS(:,:,:) = 0.0D0181:      IF(STRESST) STRESS(:,:,:) = 0.0D0
130:      !182:      !
131:      DO J1=1,NATOMS-1183:      ! Double loop over atom types
132:         J13 = 3*(J1-1)184:      LT3: DO T1=1,NSPECIES(0)
133:         T1=LABELS(J1,IP)185:         LT4: DO T2=1,T1 
134:         DO J2=J1+1,NATOMS 
135:            J23 = 3*(J2-1) 
136:            T2=LABELS(J2,IP) 
137:            ! 
138:            RHO_DUM = RHO(J1) + RHO(J2) ! factor of 1/2 is included 
139:            DUDR = REP_IJ(J1,J2) - RHO_DUM*ATT_IJ(J1,J2) 
140:            ! Want 1/R*dU/dR 
141:            DIST = RIJ(J1,J2) 
142:            DUDR = DUDR/DIST 
143:            !186:            !
144:            ! Cartesian components of distance187:            ! Double loop over free / frozen groups
145:            DO I=1,3188:            LG3: DO G1 = 1,2 ! loop over goups 1 and 2
146:               DX(I)=X(J13+I)-X(J23+I) 
147:               IGRAD(I)=DUDR*DX(I)  
148:               GRAD(J13+I) = GRAD(J13+I) + IGRAD(I) 
149:               GRAD(J23+I) = GRAD(J23+I) - IGRAD(I) 
150:            ENDDO 
151:            ! 
152:            IF(STRESST) THEN ! Accumulate local stresses 
153:               DO I=1,3 
154:                  DO J=I,3 
155:                     DUMMY = IGRAD(I)*DX(J) 
156:                     STRESS(J1,I,J) = STRESS(J1,I,J) + & 
157:                          DUMMY 
158:                     STRESS(J2,I,J) = STRESS(J2,I,J) + & 
159:                          DUMMY 
160:                     STRESS(0,I,J) = STRESS(0,I,J) + & 
161:                          DUMMY 
162:                  ENDDO 
163:               ENDDO 
164:            ENDIF 
165:            ! 
166:            IF(HESST) THEN 
167:               !189:               !
168:               DUMMY=1.0D0/DIST**2 ! Reset DUMMY190:               IF(T1 == T2) THEN
169:               REP_DUM1=REP_IJ(J1,J2)/DIST191:                  G2START = G1
170:               REP_DUM2=MGUPTA_MP_DIVBY_R0(T1,T2)*&192:               ELSE
171:                    REP_IJ(J1,J2)193:                  G2START = 1
172:               REP_DUM2=(REP_DUM2-REP_DUM1)*DUMMY194:               ENDIF
173:               ATT_DUM1=ATT_IJ(J1,J2)/DIST195:               !
174:               ATT_DUM2=MGUPTA_M2Q_DIVBY_R0(T1,T2)*&196:               LG4: DO G2 = G2START,2 ! loop over groups 1 and 2
175:                    ATT_IJ(J1,J2)197:                  !
176:               ATT_DUM2=(ATT_DUM2-ATT_DUM1)*DUMMY198:                  LI1=1
177:               !199:                  UI1=ATOMLISTS(T1,G1,0)
178:               DO I = 1,3 ! coordinates of atom J1200:                  !
179:                  DO J= 1,3 ! coordinate of atom J2 (/=J1)201:                  IF(T1==T2.AND.G1==G2) THEN
180:                     !                                202:                     SELF = .TRUE.
181:                     DUMMY = DX(I)*DX(J) ! Reset DUMMY203:                     UI1 = UI1 - 1 ! Exclude self-interaction
182:                     !204:                  ELSE
183:                     REP_DUM3 = REP_DUM2*DUMMY205:                     SELF = .FALSE.
184:                     ATT_DUM3 = ATT_DUM2*DUMMY206:                  ENDIF
185:                     !207:                  !
186:                     IF (I == J) THEN208:                  LA3: DO K1=LI1,UI1
187:                        REP_DUM3 = REP_DUM3 + REP_DUM1 
188:                        ATT_DUM3 = ATT_DUM3 + ATT_DUM1 
189:                     ENDIF 
190:                     !209:                     !
191:                     D2UDR2 = REP_DUM3 - RHO_DUM*ATT_DUM3210:                     J1 = ATOMLISTS(T1,G1,K1) ! actual atom 1 index
192:                     DUDR = D2UDR2211:                     J13 = 3*(J1-1)
193:                     !212:                     !
194:                     ! Now go through triplets... FECK!213:                     ! 2nd atom loop bounds depend on group.
195:                     DO J3=1,NATOMS                    214:                     ! Can this IF block be moved outside of loop A1?
196:                        J33 = 3*(J3-1)215:                     IF(SELF) THEN
197:                        T3 = LABELS(J3,IP)216:                        LI2=K1+1 
198:                        !217:                        UI2=UI1+1
199:                        IF ((J3.NE.J2).AND.(J3.NE.J1)) THEN218:                     ELSE
200:                           DUMMY = RHO4HESS(J3)*ATT_IJ(J1,J3)* &219:                        LI2=1
201:                                ATT_IJ(J2,J3)/(RIJ(J1,J3)*RIJ(J2,J3))220:                        UI2=ATOMLISTS(T2,G2,0)
202:                           D2UDR2 = D2UDR2 - &221:                     ENDIF
203:                                DUMMY*(X(J13+I)-X(J33+I))* &222:                     !
204:                                (X(J23+J)-X(J33+J))223:                     LA4: DO K2=LI2,UI2
205:                           DUDR = DUDR - & 
206:                                DUMMY*(X(J23+I)-X(J33+I))* & 
207:                                (X(J13+J)-X(J33+J)) 
208:                        ENDIF 
209:                        !224:                        !
210:                        IF (J3.NE.J2) THEN225:                        J2 = ATOMLISTS(T2,G2,K2) ! actual atom 2 index
211:                           DUMMY = RHO4HESS(J2)* &226:                        J23 = 3*(J2-1)
212:                                ATT_IJ(J1,J2)*ATT_IJ(J3,J2)/ & 
213:                                (RIJ(J1,J2)*RIJ(J3,J2)) 
214:                           D2UDR2 = D2UDR2 + & 
215:                                DUMMY*(X(J13+I)-X(J23+I))* & 
216:                                (X(J33+J)-X(J23+J)) 
217:                           DUDR = DUDR + & 
218:                                DUMMY*(X(J33+I)-X(J23+I))* & 
219:                                (X(J13+J)-X(J23+J)) 
220:                        ENDIF 
221:                        !227:                        !
222:                        IF (J3.NE.J1) THEN228:                        RHO_DUM = RHO(J1) + RHO(J2) ! factor of 1/2 is included
223:                           DUMMY = RHO4HESS(J1)* &229:                        DUDR = REP_IJ(J1,J2) - RHO_DUM*ATT_IJ(J1,J2)
224:                                ATT_IJ(J3,J1)*ATT_IJ(J2,J1)/ &230:                        ! Want 1/R*dU/dR
225:                                (RIJ(J3,J1)*RIJ(J2,J1))231:                        DIST = RIJ(J1,J2)
226:                           D2UDR2 = D2UDR2 + &232:                        DUDR = DUDR/DIST
227:                                DUMMY*(X(J33+I)-X(J13+I))* &233:                        !
228:                                (X(J23+J)-X(J13+J))234:                        ! Cartesian components of distance
229:                           DUDR = DUDR + &235:                        DO I=1,3
230:                                DUMMY*(X(J23+I)-X(J13+I))* &236:                           DX(I)=X(J13+I)-X(J23+I)
231:                                (X(J33+J)-X(J13+J))237:                           IGRAD(I)=DUDR*DX(I) 
 238:                           GRAD(J13+I) = GRAD(J13+I) + IGRAD(I)
 239:                           GRAD(J23+I) = GRAD(J23+I) - IGRAD(I)
 240:                        ENDDO
 241:                        !
 242:                        IF(STRESST) THEN ! Accumulate local stresses
 243:                           DO I=1,3
 244:                              DO J=I,3
 245:                                 DUMMY = IGRAD(I)*DX(J)
 246:                                 STRESS(J1,I,J) = STRESS(J1,I,J) + &
 247:                                      DUMMY
 248:                                 STRESS(J2,I,J) = STRESS(J2,I,J) + &
 249:                                      DUMMY
 250:                                 STRESS(0,I,J) = STRESS(0,I,J) + &
 251:                                      DUMMY
 252:                              ENDDO
 253:                           ENDDO
232:                        ENDIF254:                        ENDIF
233:                        !255:                        !
234:                     ENDDO256:                        IF(HESST) THEN
235:                     !257:                           !
236:                     ! Accumulate diagonal blocks first258:                           DUMMY=1.0D0/DIST**2 ! Reset DUMMY
237:                     HESS(J13+I, J13+J) = &259:                           REP_DUM1=REP_IJ(J1,J2)/DIST
238:                          HESS(J13+I, J13+J) + D2UDR2260:                           REP_DUM2=MGUPTA_MP_DIVBY_R0(T1,T2)*&
239:                     HESS(J23+I, J23+J) = &261:                                REP_IJ(J1,J2)
240:                          HESS(J23+I, J23+J) + DUDR ! D2UDR2 /= DUDR262:                           REP_DUM2=(REP_DUM2-REP_DUM1)*DUMMY
241:                     !263:                           ATT_DUM1=ATT_IJ(J1,J2)/DIST
242:                     ! Compute (not accumulate!) off-diagonal blocks,264:                           ATT_DUM2=MGUPTA_M2Q_DIVBY_R0(T1,T2)*&
243:                     ! imposing symmetry over J1 and J2.265:                                ATT_IJ(J1,J2)
244:                     D2UDR2 = -D2UDR2266:                           ATT_DUM2=(ATT_DUM2-ATT_DUM1)*DUMMY
245:                     HESS(J13+I, J23+J) = D2UDR2267:                           !
246:                     HESS(J23+J, J13+I) = D2UDR2268:                           DO I = 1,3 ! coordinates of atom J1
247:                     !269:                              DO J= 1,3 ! coordinate of atom J2 (/=J1 by construction!)
248:                  ENDDO270:                                 !                                
249:               ENDDO271:                                 DUMMY = DX(I)*DX(J) ! Reset DUMMY
250:               !272:                                 !
251:            ENDIF ! End of Hessian(J1,J2) block273:                                 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 LA4
 343:                  ENDDO LA3
 344:                  !
 345:               ENDDO LG4
 346:            ENDDO LG3
252:            !347:            !
 348:         ENDDO LT4
 349:      ENDDO LT3
 350:      !
 351:      ! Finally, zero gradiend for group 2 atoms
 352:      DO T1=1,NSPECIES(0) ! Atom types
 353:         DO J2=1,ATOMLISTS(T1,2,0) ! Group 2
 354:            J1=ATOMLISTS(T1,1,J2)  ! Actual atom index
 355:            J13=3*(J1-1)
 356:            DO I=1,3
 357:               GRAD(J13+I)=0.0D0
 358:            ENDDO
253:         ENDDO359:         ENDDO
254:      ENDDO360:      ENDDO
255:      !361:      !
256:      IF(STRESST) THEN362:      IF(STRESST) THEN
257:         DO I=1,3363:         DO I=1,3
258:            DO J=I,3364:            DO J=I,3
259:               STRESS(0,J,I) = STRESS(0,I,J) ! Impose symmetry365:               STRESS(0,J,I) = STRESS(0,I,J) ! Impose symmetry
260:            ENDDO366:            ENDDO
261:         ENDDO367:         ENDDO
262:      ENDIF368:      ENDIF


r31119/MLJ.f90 2016-09-06 15:30:36.922417925 +0100 r31118/MLJ.f90 2016-09-06 15:30:40.142460898 +0100
 19: !  =============================================================== 19: !  ===============================================================
 20: !  Multicomponent Lennard-Jones without a cutoff. 20: !  Multicomponent Lennard-Jones without a cutoff.
 21: !  Assumed reduced units with SIGMA_AA = EPSILON_AA = 1. 21: !  Assumed reduced units with SIGMA_AA = EPSILON_AA = 1.
 22: !  Per-atom energy is stored in array VT(NATOMS). 22: !  Per-atom energy is stored in array VT(NATOMS).
 23: ! 23: !
 24: !  ds656> 6/12/2014 24: !  ds656> 6/12/2014
 25: !  =============================================================== 25: !  ===============================================================
 26: ! 26: !
 27: SUBROUTINE MLJ(X,GRAD,POT,GRADT,HESST,STRESST) 27: SUBROUTINE MLJ(X,GRAD,POT,GRADT,HESST,STRESST)
 28:   ! 28:   !
 29:   USE COMMONS, ONLY : NATOMS,VT,LABELS,NSPECIES,STRESS,MYNODE 29:   USE COMMONS, ONLY : NATOMS,VT,ATOMLISTS, NSPECIES, STRESS
 30:   USE MODHESS, ONLY : HESS 30:   USE MODHESS, ONLY : HESS
 31:   USE POT_PARAMS, ONLY : MLJ_EPS, MLJ_SIG 31:   USE POT_PARAMS, ONLY : MLJ_EPS, MLJ_SIG
 32:   ! 32:   !
 33:   IMPLICIT NONE 33:   IMPLICIT NONE
 34:   ! 34:   !
 35:   LOGICAL, INTENT (IN) :: GRADT,HESST,STRESST 35:   LOGICAL, INTENT (IN) :: GRADT,HESST,STRESST
 36:   DOUBLE PRECISION, INTENT (IN) :: X(3*NATOMS) 36:   DOUBLE PRECISION, INTENT (IN) :: X(3*NATOMS)
 37:   DOUBLE PRECISION, INTENT (INOUT) :: POT, GRAD(3*NATOMS) 37:   DOUBLE PRECISION, INTENT (INOUT) :: POT, GRAD(3*NATOMS)
 38:   ! 38:   !
 39:   INTEGER :: I,J,J1,J13,J2,J23,T1,T2,IP 39:   LOGICAL :: SELF
  40:   INTEGER :: I,J,J1,J13,J2,J23,T1,T2,G1,G2,G2START,&
  41:        K1,K2,LI1,LI2,UI1,UI2
 40:   DOUBLE PRECISION :: DIST2, DX(3), IR6, IR12, DUMMY, & 42:   DOUBLE PRECISION :: DIST2, DX(3), IR6, IR12, DUMMY, &
 41:        SIG2(NSPECIES(0),NSPECIES(0)), & 43:        SIG2, EPS4, EPS24, IGRAD(3), DUDR, D2UDR2
 42:        EPS4(NSPECIES(0),NSPECIES(0)), & 
 43:        EPS24(NSPECIES(0),NSPECIES(0)), IGRAD(3), DUDR, D2UDR2 
 44:   ! 44:   !
 45:   ! Zero the potential and the gradient 45:   ! Zero the potential and the gradient
 46:   VT(1:NATOMS) = 0.0D0 46:   VT(1:NATOMS) = 0.0D0
 47:   POT = 0.0D0 47:   POT = 0.0D0
 48:   IF(GRADT) GRAD(:) = 0.0D0 48:   IF(GRADT) GRAD(:) = 0.0D0
 49:   IF(HESST) HESS(:,:) = 0.0D0 49:   IF(HESST) HESS(:,:) = 0.0D0
 50:   IF(STRESST) STRESS(:,:,:) = 0.0D0 50:   IF(STRESST) STRESS(:,:,:) = 0.0D0
 51:   ! 51:   !
 52:   ! This can be precomputed elsewhere... 
 53:   SIG2(:,:) = MLJ_SIG(:,:)*MLJ_SIG(:,:) 
 54:   EPS4(:,:) = 4.0D0*MLJ_EPS(:,:) 
 55:   EPS24(:,:) = 24.0D0*MLJ_EPS(:,:) 
 56:   ! 
 57:   IP=MYNODE+1 
 58:   ! 
 59:   ! Double loop over atom types 52:   ! Double loop over atom types
 60:   DO J1=1,NATOMS-1 53:   LT1: DO T1=1,NSPECIES(0)
 61:      J13=3*(J1-1) 54:      LT2: DO T2=1,T1 
 62:      T1=LABELS(J1,IP) 
 63:      DO J2=J1+1,NATOMS 
 64:         J23 = 3*(J2-1) 
 65:         T2=LABELS(J2,IP) 
 66:         ! 
 67:         ! 
 68:         DIST2 = 0.0D0 
 69:         DO I = 1,3 
 70:            DX(I) = X(J13 + I) - X(J23 + I) 
 71:            DIST2 = DIST2 + DX(I)*DX(I) 
 72:         ENDDO 
 73:         ! 55:         !
 74:         ! ---- Potential-specific bit ------------- 56:         ! This can be precomputed elsewhere...
 75:         IR6 = (SIG2(T1,T2)/DIST2)**3 57:         SIG2 = MLJ_SIG(T1,T2)*MLJ_SIG(T1,T2)
 76:         IR12 = IR6*IR6 58:         EPS4 = 4.0D0*MLJ_EPS(T1,T2)
 77:         DUMMY = EPS4(T1,T2)*(IR12 - IR6) 59:         IF(GRADT) EPS24 = 24.0D0*MLJ_EPS(T1,T2)
 78:         !WRITE(*,*) "DUMMY FOR POT", DUMMY 
 79:         ! ----------------------------------------- 
 80:         ! 60:         !
 81:         VT(J1) = VT(J1) + DUMMY 61:         ! Double loop over mobile and frozen atom groups
 82:         VT(J2) = VT(J2) + DUMMY 62:         LG1: DO G1 = 1,2 ! 1 -> mobile; 2 -> frozen
 83:         POT = POT + DUMMY 
 84:         ! 
 85:         IF(GRADT) THEN ! Calculate gradient 
 86:            ! Compute 1/R*dUdR 
 87:            DUDR = EPS24(T1,T2)*(IR6 - 2.0D0*IR12)/DIST2 
 88:            !WRITE(*,*) "DUMMY FOR GRAD", DUMMY 
 89:            DO I = 1,3 
 90:               IGRAD(I) = DUDR*DX(I) 
 91:               GRAD(J13+I) = GRAD(J13+I) + IGRAD(I) 
 92:               GRAD(J23+I) = GRAD(J23+I) - IGRAD(I) 
 93:            ENDDO 
 94:            ! 63:            !
 95:            IF(STRESST) THEN 64:            IF(T1 == T2) THEN
 96:               DO I=1,3 65:               G2START = G1
 97:                  DO J=I,3 66:            ELSE
 98:                     DUMMY = IGRAD(I)*DX(J) 67:               G2START = 1
 99:                     STRESS(J1,I,J) = STRESS(J1,I,J) + & 
100:                          DUMMY 
101:                     STRESS(J2,I,J) = STRESS(J2,I,J) + & 
102:                          DUMMY 
103:                     STRESS(0,I,J) = STRESS(0,I,J) + & 
104:                          DUMMY 
105:                  ENDDO 
106:               ENDDO 
107:            ENDIF 68:            ENDIF
108:            ! 69:            !
109:            IF(HESST) THEN 70:            LG2: DO G2 = G2START,2 ! include frozen-frozen
110:               ! Compute d2UdR2 71:               !
111:               D2UDR2 = EPS24(T1,T2)*(26.0D0*IR12-7.0D0*IR6)/DIST2 72:               LI1=1
112:               ! What we need is (D2UDR - DUDR/DIST) / DIST 73:               UI1=ATOMLISTS(T1,G1,0)
113:               D2UDR2 = (D2UDR2 - DUDR) / DIST2 
114:               ! 74:               !
115:               DO I=1,3 75:               IF(T1==T2.AND.G1==G2) THEN
116:                  DO J=I,3 76:                  SELF = .TRUE.
  77:                  UI1 = UI1 - 1 ! exclude self-interaction
  78:               ELSE
  79:                  SELF = .FALSE.
  80:               ENDIF
  81:               !
  82:               LA1: DO K1=LI1,UI1
  83:                  !
  84:                  J1 = ATOMLISTS(T1,G1,K1) ! actual atom 1 index
  85:                  J13 = 3*(J1-1)
  86:                  !
  87:                  ! Atom loop bounds depend on type and group.
  88:                  ! Can this IF block be moved outside of loop LA1?
  89:                  IF(SELF) THEN
  90:                     LI2=K1+1
  91:                     UI2=UI1+1
  92:                  ELSE
  93:                     LI2=1
  94:                     UI2=ATOMLISTS(T2,G2,0)
  95:                  ENDIF
  96:                  !
  97:                  LA2: DO K2=LI2,UI2
  98:                     !
  99:                     J2 = ATOMLISTS(T2,G2,K2) ! actual atom 2 index
 100:                     J23 = 3*(J2-1)
 101:                     !
 102:                     DIST2 = 0.0D0
 103:                     DO I = 1,3
 104:                        DX(I) = X(J13 + I) - X(J23 + I)
 105:                        DIST2 = DIST2 + DX(I)*DX(I)
 106:                     ENDDO
 107:                     !
 108:                     ! ---- Potential-specific bit -------------
 109:                     IR6 = (SIG2/DIST2)**3
 110:                     IR12 = IR6*IR6
 111:                     DUMMY = EPS4*(IR12 - IR6)
 112:                     !WRITE(*,*) "DUMMY FOR POT", DUMMY
 113:                     ! -----------------------------------------
117:                     !114:                     !
118:                     DUMMY = DX(I)*DX(J) ! Reset DUMMY115:                     VT(J1) = VT(J1) + DUMMY
119:                     DUMMY = DUMMY*D2UDR2116:                     VT(J2) = VT(J2) + DUMMY
120:                     IF(I==J) DUMMY = DUMMY + DUDR117:                     POT = POT + DUMMY
121:                     !118:                     !
122:                     ! Accumulate diagonal blocks first119:                     IF(GRADT) THEN ! Calculate gradient
123:                     HESS(J13+I, J13+J) = &120:                        ! Compute 1/R*dUdR
124:                          HESS(J13+I, J13+J) + DUMMY121:                        DUDR = EPS24*(IR6 - 2.0D0*IR12)/DIST2
125:                     HESS(J23+I, J23+J) = &122:                        !WRITE(*,*) "DUMMY FOR GRAD", DUMMY
126:                          HESS(J23+I, J23+J) + DUMMY123:                        DO I = 1,3
127:                     ! NOTE: J < I will be filled later 124:                           IGRAD(I) = DUDR*DX(I)
128:                     ! by symmetry of Hessian125:                           GRAD(J13+I) = GRAD(J13+I) + IGRAD(I)
129:                     !                                126:                           GRAD(J23+I) = GRAD(J23+I) - IGRAD(I)
130:                     ! Off-diagonal blocks acquire 127:                        ENDDO
131:                     ! contribution from only one128:                        !
132:                     ! pair (J1,J2) in pairwise-additive 129:                        IF(STRESST) THEN
133:                     ! potentials, so no  130:                           DO I=1,3
134:                     ! accumulation is needed 131:                              DO J=I,3
135:                     ! (just value assignment).132:                                 DUMMY = IGRAD(I)*DX(J)
136:                     !133:                                 STRESS(J1,I,J) = STRESS(J1,I,J) + &
137:                     DUMMY = -DUMMY134:                                      DUMMY
138:                     HESS(J13+I, J23+J) = DUMMY135:                                 STRESS(J2,I,J) = STRESS(J2,I,J) + &
139:                     HESS(J23+J, J13+I) = DUMMY 136:                                      DUMMY
140:                     ! imposed symmetry 137:                                 STRESS(0,I,J) = STRESS(0,I,J) + &
141:                     !138:                                      DUMMY
142:                     IF(I .NE. J) THEN ! Use symmetry of DUMMY139:                              ENDDO
143:                        HESS(J13+J, J23+I) = DUMMY140:                           ENDDO
144:                        HESS(J23+I, J13+J) = DUMMY141:                        ENDIF
 142:                        !
 143:                        IF(HESST) THEN
 144:                           ! Compute d2UdR2
 145:                           D2UDR2 = EPS24*(26.0D0*IR12-7.0D0*IR6)/DIST2
 146:                           ! What we need is (D2UDR - DUDR/DIST) / DIST
 147:                           D2UDR2 = (D2UDR2 - DUDR) / DIST2
 148:                           !
 149:                           DO I=1,3
 150:                              DO J=I,3
 151:                                 !
 152:                                 DUMMY = DX(I)*DX(J) ! Reset DUMMY
 153:                                 DUMMY = DUMMY*D2UDR2
 154:                                 IF(I==J) DUMMY = DUMMY + DUDR
 155:                                 !
 156:                                 ! Accumulate diagonal blocks first
 157:                                 HESS(J13+I, J13+J) = &
 158:                                      HESS(J13+I, J13+J) + DUMMY
 159:                                 HESS(J23+I, J23+J) = &
 160:                                      HESS(J23+I, J23+J) + DUMMY
 161:                                 ! NOTE: J < I will be filled later 
 162:                                 ! by symmetry of Hessian
 163:                                 !                                
 164:                                 ! Off-diagonal blocks acquire 
 165:                                 ! contribution from only one
 166:                                 ! pair (J1,J2) in pairwise-additive 
 167:                                 ! potentials, so no  
 168:                                 ! accumulation is needed 
 169:                                 ! (just value assignment).
 170:                                 !
 171:                                 DUMMY = -DUMMY
 172:                                 HESS(J13+I, J23+J) = DUMMY
 173:                                 HESS(J23+J, J13+I) = DUMMY 
 174:                                 ! imposed symmetry 
 175:                                 !
 176:                                 IF(I .NE. J) THEN ! Use symmetry of DUMMY w.r.t. K,L.
 177:                                    HESS(J13+J, J23+I) = DUMMY
 178:                                    HESS(J23+I, J13+J) = DUMMY
 179:                                 ENDIF
 180:                                 !
 181:                              ENDDO
 182:                           ENDDO
 183:                           !
 184:                        ENDIF ! End of HESSIAN loop
 185:                        !
145:                     ENDIF186:                     ENDIF
146:                     !187:                     !
147:                  ENDDO188:                  ENDDO LA2 ! Double loop over atoms
148:               ENDDO189:               ENDDO LA1
149:               !190:               !
150:            ENDIF ! End of HESSIAN loop191:            ENDDO LG2 ! Double loop over groups
151:            !192:         ENDDO LG1
152:         ENDIF 
153:         !193:         !
154:      ENDDO  ! Double loop over atoms194:      ENDDO LT2 ! Double loop over types
155:   ENDDO195:   ENDDO LT1
156:   !196:   !
157:   IF(GRADT) THEN197:   IF(GRADT) THEN
 198:      DO T1=1,NSPECIES(0) ! Atom types
 199:         DO J2=1,ATOMLISTS(T1,2,0) ! span group 2 only
 200:            J1=ATOMLISTS(T1,2,J2)
 201:            J13=3*(J1-1)
 202:            GRAD(J13+1) = 0.0D0 ! Reset to zero
 203:            GRAD(J13+2) = 0.0D0
 204:            GRAD(J13+3) = 0.0D0
 205:         ENDDO
 206:      ENDDO
158:      !207:      !
159:      ! Test per-atom energies and gradients:208:      ! Test per-atom energies and gradients:
160:      ! DO T1=1,NSPECIES(0)209:      ! DO T1=1,NSPECIES(0)
161:      !    DO G1=1,2210:      !    DO G1=1,2
162:      !       WRITE(*,'(A,3(1X,I5))') "BLJ_CLUST> T,G,NAT(T,G)=", &211:      !       WRITE(*,'(A,3(1X,I5))') "BLJ_CLUST> T,G,NAT(T,G)=", &
163:      !            T1, G1, ATOMLISTS(T1,G1,0)212:      !            T1, G1, ATOMLISTS(T1,G1,0)
164:      !       DO J2=1,ATOMLISTS(T1,G1,0)213:      !       DO J2=1,ATOMLISTS(T1,G1,0)
165:      !          J1=ATOMLISTS(T1,G1,J2)214:      !          J1=ATOMLISTS(T1,G1,J2)
166:      !          J13=3*(J1-1)215:      !          J13=3*(J1-1)
167:      !          WRITE(*,'(A,1X, I4, 4(1X,F12.6))') "BLJ_CLUST>", &216:      !          WRITE(*,'(A,1X, I4, 4(1X,F12.6))') "BLJ_CLUST>", &
172:      !221:      !
173:      IF(STRESST) THEN222:      IF(STRESST) THEN
174:         DO I=1,3223:         DO I=1,3
175:            DO J=I,3224:            DO J=I,3
176:               STRESS(0,J,I) = STRESS(0,I,J) ! Impose symmetry225:               STRESS(0,J,I) = STRESS(0,I,J) ! Impose symmetry
177:            ENDDO226:            ENDDO
178:         ENDDO227:         ENDDO
179:      ENDIF228:      ENDIF
180:      !229:      !
181:      IF(HESST) THEN230:      IF(HESST) THEN
182:         ! Use symmetry to fill skipped entries in diagonal blocks231:         ! Use symmetry to fill the skipped entries in diagonal blocks
183:         DO J1=1,NATOMS232:         DO J1=1,NATOMS
184:            DO I = 1,3 ! coordinates of atom 1233:            DO I = 1,3 ! coordinates of atom 1
185:               DO J = I+1,3 ! coordinates of atom 2 234:               DO J = I+1,3 ! coordinates of atom 2 
186:                  HESS(J1+J, J1+I) = HESS(J1+I, J1+J)235:                  HESS(J1+J, J1+I) = HESS(J1+I, J1+J)
187:               ENDDO236:               ENDDO
188:            ENDDO237:            ENDDO
189:         ENDDO238:         ENDDO
190:         !239:         !
191:      ENDIF240:      ENDIF
192:      !241:      !


r31119/QALCSearch.f90 2016-09-06 15:30:37.190421496 +0100 r31118/QALCSearch.f90 2016-09-06 15:30:40.442464900 +0100
376:      DO J=1, LISTS_AB(2,0)376:      DO J=1, LISTS_AB(2,0)
377:         !377:         !
378:         ! Append to list of swaps only if swap is non-degenerate378:         ! Append to list of swaps only if swap is non-degenerate
379:         !379:         !
380:         N = N + 1380:         N = N + 1
381:         SWAPS(N,1) = LISTS_AB(1,I)381:         SWAPS(N,1) = LISTS_AB(1,I)
382:         SWAPS(N,2) = LISTS_AB(2,J)382:         SWAPS(N,2) = LISTS_AB(2,J)
383:         !383:         !
384:         IF(QALCSMODE==3) THEN ! Sort list384:         IF(QALCSMODE==3) THEN ! Sort list
385:            ! Perform swap and compute unquenched energy.385:            ! Perform swap and compute unquenched energy.
386:            CALL SWAP_LABELS(LISTS_AB(1,I),LISTS_AB(2,J),NP)386:            CALL SWAP_LABELS(LISTS_AB(1,I),LISTS_AB(2,J))
387:            CALL POTENTIAL(COORDS(:,NP),X0(:),EST(N),.FALSE.,.FALSE.)387:            CALL POTENTIAL(COORDS(:,NP),X0(:),EST(N),.FALSE.,.FALSE.)
388:            CALL SWAP_LABELS(LISTS_AB(1,I),LISTS_AB(2,J),NP)388:            CALL SWAP_LABELS(LISTS_AB(1,I),LISTS_AB(2,J))
389:            ! Sort by unquenched energy389:            ! Sort by unquenched energy
390:            SORT_EST: DO K=N,2,-1390:            SORT_EST: DO K=N,2,-1
391:               IF(EST(K) < EST(K-1)) THEN391:               IF(EST(K) < EST(K-1)) THEN
392:                  E0 = EST(K); EST(K) = EST(K-1); EST(K-1) = E0392:                  E0 = EST(K); EST(K) = EST(K-1); EST(K-1) = E0
393:                  I_AB(1:2) = SWAPS(K,1:2)393:                  I_AB(1:2) = SWAPS(K,1:2)
394:                  SWAPS(K,1:2) = SWAPS(K-1,1:2)394:                  SWAPS(K,1:2) = SWAPS(K-1,1:2)
395:                  SWAPS(K-1,1:2) = I_AB(1:2)395:                  SWAPS(K-1,1:2) = I_AB(1:2)
396:               ELSE396:               ELSE
397:                  EXIT SORT_EST397:                  EXIT SORT_EST
398:               ENDIF398:               ENDIF
411:   !411:   !
412:   NTOT = N412:   NTOT = N
413:   DO I=1,NTOT413:   DO I=1,NTOT
414:      !414:      !
415:      IF(QALCSMODE==2) THEN    415:      IF(QALCSMODE==2) THEN    
416:         J = INT(DPRAND()*DBLE(N)) + 1416:         J = INT(DPRAND()*DBLE(N)) + 1
417:      ELSEIF(QALCSMODE==3) THEN417:      ELSEIF(QALCSMODE==3) THEN
418:         J=I418:         J=I
419:      ENDIF419:      ENDIF
420:      !420:      !
421:      CALL SWAP_LABELS(SWAPS(J,1),SWAPS(J,2),NP)421:      CALL SWAP_LABELS(SWAPS(J,1),SWAPS(J,2))
422:      !                                                            422:      !                                                            
423:      NQTOT = NQTOT + 1423:      NQTOT = NQTOT + 1
424:      NQ(NP) = NQ(NP) + 1424:      NQ(NP) = NQ(NP) + 1
425:      CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)425:      CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)
426:      IF(QALCSV) THEN426:      IF(QALCSV) THEN
427:         CALL PRINT_QUENCH(NP, ITER, '  ')427:         CALL PRINT_QUENCH(NP, ITER, '  ')
428:         IF(QALCSMODE==3) THEN428:         IF(QALCSMODE==3) THEN
429:            WRITE(MYUNIT,'(2(A,F15.10))') &429:            WRITE(MYUNIT,'(2(A,F15.10))') &
430:                 'ab_search2> dE*= ',EST(J)-E0,' dE= ',POTEL-E0 430:                 'ab_search2> dE*= ',EST(J)-E0,' dE= ',POTEL-E0 
431:         ENDIF431:         ENDIF
433:      !433:      !
434:      IF(POTEL < E0 - ECONV) THEN434:      IF(POTEL < E0 - ECONV) THEN
435:         SEQLENGTH = SEQLENGTH + 1435:         SEQLENGTH = SEQLENGTH + 1
436:         WRITE(MYUNIT,'(A,I4,A,I4,A,F15.10)') &436:         WRITE(MYUNIT,'(A,I4,A,I4,A,F15.10)') &
437:              'ab_search2> ',SWAPS(J,1),' <-> ',SWAPS(J,2),&437:              'ab_search2> ',SWAPS(J,1),' <-> ',SWAPS(J,2),&
438:              ' => dE= ', POTEL-E0438:              ' => dE= ', POTEL-E0
439:         I_AB(1:2) = SWAPS(J,1:2)439:         I_AB(1:2) = SWAPS(J,1:2)
440:         RETURN440:         RETURN
441:      ELSE 441:      ELSE 
442:         ! Undo the swap and restore the configuratiion442:         ! Undo the swap and restore the configuratiion
443:         CALL SWAP_LABELS(SWAPS(J,1),SWAPS(J,2),NP)443:         CALL SWAP_LABELS(SWAPS(J,1),SWAPS(J,2))
444:         COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)444:         COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)
445:         POTEL = E0445:         POTEL = E0
446:         !446:         !
447:         IF(QALCSMODE==2) THEN447:         IF(QALCSMODE==2) THEN
448:            ! Remove the rejected swap from list448:            ! Remove the rejected swap from list
449:            SWAPS(J,1:2) = SWAPS(N,1:2)449:            SWAPS(J,1:2) = SWAPS(N,1:2)
450:            SWAPS(N,1:2) = 0450:            SWAPS(N,1:2) = 0
451:            N = N - 1451:            N = N - 1
452:         ENDIF452:         ENDIF
453:         !453:         !
558:         !<ds656 ...testing558:         !<ds656 ...testing
559:         !559:         !
560:         scan: DO IPICK1=1,SUBLISTS_AB(TPICK1,0)560:         scan: DO IPICK1=1,SUBLISTS_AB(TPICK1,0)
561:            ! 561:            ! 
562:            !WRITE(MYUNIT,*) 'AB_SWAP_SEARCH> TPICK1,IPICK1,I,:', &562:            !WRITE(MYUNIT,*) 'AB_SWAP_SEARCH> TPICK1,IPICK1,I,:', &
563:            !     TPICK1, IPICK1, SUBLISTS_AB(TPICK1, IPICK1)563:            !     TPICK1, IPICK1, SUBLISTS_AB(TPICK1, IPICK1)
564:            !WRITE(MYUNIT,*) 'AB_SWAP_SEARCH> TPICK2,IPICK2,J,:', &564:            !WRITE(MYUNIT,*) 'AB_SWAP_SEARCH> TPICK2,IPICK2,J,:', &
565:            !     TPICK2, IPICK2, SUBLISTS_AB(TPICK2, IPICK2)565:            !     TPICK2, IPICK2, SUBLISTS_AB(TPICK2, IPICK2)
566:            !566:            !
567:            CALL SWAP_LABELS( SUBLISTS_AB(TPICK1,IPICK1), &567:            CALL SWAP_LABELS( SUBLISTS_AB(TPICK1,IPICK1), &
568:                 SUBLISTS_AB(TPICK2,IPICK2), NP )568:                 SUBLISTS_AB(TPICK2,IPICK2) )
569:            !569:            !
570:            NQTOT = NQTOT + 1570:            NQTOT = NQTOT + 1
571:            NQ(NP) = NQ(NP) + 1571:            NQ(NP) = NQ(NP) + 1
572:            CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)572:            CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)
573:            IF(QALCSV) CALL PRINT_QUENCH(NP, ITER, '  ')573:            IF(QALCSV) CALL PRINT_QUENCH(NP, ITER, '  ')
574:            !574:            !
575:            IF(POTEL < EBEST - ECONV) THEN575:            IF(POTEL < EBEST - ECONV) THEN
576:               COMPLETE = .TRUE.576:               COMPLETE = .TRUE.
577:               STUCK = .FALSE.577:               STUCK = .FALSE.
578:               ! Store best-encountered swap candidates578:               ! Store best-encountered swap candidates
582:               EBEST = POTEL582:               EBEST = POTEL
583:               XBEST(1:3*NATOMS) = COORDS(1:3*NATOMS,NP)583:               XBEST(1:3*NATOMS) = COORDS(1:3*NATOMS,NP)
584:               WRITE(MYUNIT,'(A,I4,A,I4,A,F15.10)') &584:               WRITE(MYUNIT,'(A,I4,A,I4,A,F15.10)') &
585:                    'ab_search> ', I_AB(1),' <-> ', I_AB(2), &585:                    'ab_search> ', I_AB(1),' <-> ', I_AB(2), &
586:                    ' => dE= ',EBEST-E0586:                    ' => dE= ',EBEST-E0
587:               CALL FLUSH(MYUNIT)587:               CALL FLUSH(MYUNIT)
588:            ENDIF588:            ENDIF
589:            !589:            !
590:            ! Repeat the swap (to undo it) and reset the coordinates590:            ! Repeat the swap (to undo it) and reset the coordinates
591:            CALL SWAP_LABELS( SUBLISTS_AB(TPICK1,IPICK1), &591:            CALL SWAP_LABELS( SUBLISTS_AB(TPICK1,IPICK1), &
592:                 SUBLISTS_AB(TPICK2,IPICK2), NP )592:                 SUBLISTS_AB(TPICK2,IPICK2) )
593:            !SCREENC(1:3*NATOMS) = X0(1:3*NATOMS) !?593:            !SCREENC(1:3*NATOMS) = X0(1:3*NATOMS) !?
594:            COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)594:            COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)
595:            POTEL = E0595:            POTEL = E0
596:            !596:            !
597:         ENDDO scan597:         ENDDO scan
598:         !WRITE(MYUNIT,'(A)') 'ab_search> ...finished row/column!'598:         !WRITE(MYUNIT,'(A)') 'ab_search> ...finished row/column!'
599:         !599:         !
600:         IPICK1 = IPICK1BEST600:         IPICK1 = IPICK1BEST
601:         !601:         !
602:         ! Remove TPICK2, IPICK2 from SUBLISTS_AB and SUBFLIPS_AB602:         ! Remove TPICK2, IPICK2 from SUBLISTS_AB and SUBFLIPS_AB
615:         ENDIF615:         ENDIF
616:         !616:         !
617:      ENDDO sublist617:      ENDDO sublist
618:      !618:      !
619:      IF(.NOT. SCAN_ALL) COMPLETE = .TRUE.619:      IF(.NOT. SCAN_ALL) COMPLETE = .TRUE.
620:      !620:      !
621:   ENDDO list621:   ENDDO list
622:   !622:   !
623:   IF(EBEST < E0 - ECONV) THEN623:   IF(EBEST < E0 - ECONV) THEN
624:      SEQLENGTH = SEQLENGTH + 1624:      SEQLENGTH = SEQLENGTH + 1
625:      CALL SWAP_LABELS( I_AB(1), I_AB(2), NP )625:      CALL SWAP_LABELS( I_AB(1), I_AB(2) )
626:      !SCREENC(1:3*NATOMS) = XBEST(1:3*NATOMS)626:      !SCREENC(1:3*NATOMS) = XBEST(1:3*NATOMS)
627:      COORDS(1:3*NATOMS, NP) = XBEST(1:3*NATOMS)627:      COORDS(1:3*NATOMS, NP) = XBEST(1:3*NATOMS)
628:      POTEL = EBEST628:      POTEL = EBEST
629:      WRITE(MYUNIT,'(A,I4,A,I4,A,F20.10)') &629:      WRITE(MYUNIT,'(A,I4,A,I4,A,F20.10)') &
630:           'ab_search> Swapped ',I_AB(1),' and ',I_AB(2), &630:           'ab_search> Swapped ',I_AB(1),' and ',I_AB(2), &
631:           ' to get E=',POTEL631:           ' to get E=',POTEL
632:      CALL FLUSH(MYUNIT)632:      CALL FLUSH(MYUNIT)
633:   !ELSE633:   !ELSE
634:   !   WRITE(MYUNIT,'(A)') 'ab_search> Found no good swaps.'634:   !   WRITE(MYUNIT,'(A)') 'ab_search> Found no good swaps.'
635:   ENDIF635:   ENDIF
668:   !668:   !
669:   FLIPS_AB(1:2,1:NATOMS) = 0.0D0669:   FLIPS_AB(1:2,1:NATOMS) = 0.0D0
670:   X0(1:3*NATOMS) = COORDS(1:3*NATOMS, NP)670:   X0(1:3*NATOMS) = COORDS(1:3*NATOMS, NP)
671:   E0 = POTEL671:   E0 = POTEL
672:   !672:   !
673:   DO T=1,2673:   DO T=1,2
674:      DO I=1,LISTS_AB(T,0)674:      DO I=1,LISTS_AB(T,0)
675:         !675:         !
676:         ! Flip the label of current atom index and quench676:         ! Flip the label of current atom index and quench
677:         J=LISTS_AB(T,I) ! permanent atom index677:         J=LISTS_AB(T,I) ! permanent atom index
678:         CALL FLIP_LABEL(J,BA(T),NP)678:         CALL FLIP_LABEL(J,BA(T))
679:         !write(MYUNIT,'(A,4(1X,I2))') 'calc_flips> T,I,J,BA(T)=',T,I,J,BA(T)679:         !write(MYUNIT,'(A,4(1X,I2))') 'calc_flips> T,I,J,BA(T)=',T,I,J,BA(T)
680:         !680:         !
681:         TEMPSAVE = SAVEMULTIMINONLY681:         TEMPSAVE = SAVEMULTIMINONLY
682:         SAVEMULTIMINONLY = .TRUE. ! never save minima with 1 atom flipped!682:         SAVEMULTIMINONLY = .TRUE. ! never save minima with 1 atom flipped!
683:         NQTOT = NQTOT + 1683:         NQTOT = NQTOT + 1
684:         NQ(NP) = NQ(NP) + 1684:         NQ(NP) = NQ(NP) + 1
685:         CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)685:         CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)
686:         IF(QALCSV) CALL PRINT_QUENCH(NP, ITER, '* ')686:         IF(QALCSV) CALL PRINT_QUENCH(NP, ITER, '* ')
687:         SAVEMULTIMINONLY = TEMPSAVE687:         SAVEMULTIMINONLY = TEMPSAVE
688:         !688:         !
689:         ! Store the new energy after quenched flip689:         ! Store the new energy after quenched flip
690:         FLIPS_AB(T,I) = POTEL690:         FLIPS_AB(T,I) = POTEL
691:         !691:         !
692:         ! Undo the flip and revert the configuration692:         ! Undo the flip and revert the configuration
693:         CALL FLIP_LABEL(J,AB(T), NP)693:         CALL FLIP_LABEL(J,AB(T))
694:         COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)694:         COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)
695:         POTEL = E0695:         POTEL = E0
696:         !696:         !
697:      ENDDO697:      ENDDO
698:   ENDDO698:   ENDDO
699:   !699:   !
700:   !STOP ! TESTING!700:   !STOP ! TESTING!
701:   !701:   !
702:   RETURN702:   RETURN
703:   !703:   !
704: END SUBROUTINE CALC_FLIPS704: END SUBROUTINE CALC_FLIPS
705: !705: !
706: !=============================================================706: !=============================================================
707: !707: !
708: SUBROUTINE SWAP_LABELS(I,J,NP)708: SUBROUTINE SWAP_LABELS(I,J)
709:   !709:   !
710:   ! Swap the labels/types of two particles identified by their710:   ! Swap the labels/types of two particles identified by their
711:   ! permanent indices I and J.711:   ! permanent indices I and J.
712:   !712:   !
713:   ! USE COMMONS, ONLY : ATOMLISTS, INVATOMLISTS, MYUNIT, ATMASS, READMASST713:   ! USE COMMONS, ONLY : ATOMLISTS, INVATOMLISTS, MYUNIT, ATMASS, READMASST
714:   USE COMMONS, ONLY : ATOMLISTS, INVATOMLISTS, ATMASS, &714:   USE COMMONS, ONLY : ATOMLISTS, INVATOMLISTS, ATMASS, &
715:        READMASST, SPECMASST, LABELS715:        READMASST, SPECMASST
716:   !716:   !
717:   IMPLICIT NONE717:   IMPLICIT NONE
718:   !718:   !
719:   INTEGER, INTENT(IN) :: I,J,NP719:   INTEGER, INTENT(IN) :: I,J
720:   INTEGER :: ITYPE,JTYPE,IPLACE,JPLACE,IGROUP,JGROUP,K720:   INTEGER :: ITYPE,JTYPE,IPLACE,JPLACE,IGROUP,JGROUP,K
721:   DOUBLE PRECISION :: DUMMY721:   DOUBLE PRECISION :: DUMMY
722:   !722:   !
723:   !WRITE(MYUNIT,*) 'SWAPPING I,J=', I,J723:   !WRITE(MYUNIT,*) 'SWAPPING I,J=', I,J
724:   !724:   !
725:   ITYPE = INVATOMLISTS(I,1)725:   ITYPE = INVATOMLISTS(I,1)
726:   IGROUP = INVATOMLISTS(I,2)726:   IGROUP = INVATOMLISTS(I,2)
727:   IPLACE = INVATOMLISTS(I,3)727:   IPLACE = INVATOMLISTS(I,3)
728:   !728:   !
729:   JTYPE = INVATOMLISTS(J,1)729:   JTYPE = INVATOMLISTS(J,1)
730:   JGROUP = INVATOMLISTS(J,2)730:   JGROUP = INVATOMLISTS(J,2)
731:   JPLACE = INVATOMLISTS(J,3)731:   JPLACE = INVATOMLISTS(J,3)
732:   !732:   !
733:   LABELS(I,NP) = JTYPE 
734:   LABELS(J,NP) = ITYPE 
735:   ! 
736:   !WRITE(MYUNIT,*) 'ITYPE, IGROUP, IPLACE=', ITYPE, IGROUP, IPLACE733:   !WRITE(MYUNIT,*) 'ITYPE, IGROUP, IPLACE=', ITYPE, IGROUP, IPLACE
737:   !WRITE(MYUNIT,*) 'JTYPE, JGROUP, JPLACE=', JTYPE, JGROUP, JPLACE734:   !WRITE(MYUNIT,*) 'JTYPE, JGROUP, JPLACE=', JTYPE, JGROUP, JPLACE
738:   !735:   !
739:   !WRITE(MYUNIT,*) 'INITIAL ITYPE ( ', ITYPE, ' ) LIST:', &736:   !WRITE(MYUNIT,*) 'INITIAL ITYPE ( ', ITYPE, ' ) LIST:', &
740:   !     (ATOMLISTS(ITYPE,IGROUP,K), K=1,ATOMLISTS(ITYPE,IGROUP,0))737:   !     (ATOMLISTS(ITYPE,IGROUP,K), K=1,ATOMLISTS(ITYPE,IGROUP,0))
741:   !WRITE(MYUNIT,*) 'INITIAL JTYPE ( ', JTYPE, ' ) LIST:', &738:   !WRITE(MYUNIT,*) 'INITIAL JTYPE ( ', JTYPE, ' ) LIST:', &
742:   !     (ATOMLISTS(JTYPE,JGROUP,K), K=1,ATOMLISTS(JTYPE,JGROUP,0))739:   !     (ATOMLISTS(JTYPE,JGROUP,K), K=1,ATOMLISTS(JTYPE,JGROUP,0))
743:   !740:   !
744:   ! Swap labels741:   ! Swap labels
745:   ATOMLISTS(ITYPE, IGROUP, IPLACE) = J742:   ATOMLISTS(ITYPE, IGROUP, IPLACE) = J
762:   !WRITE(MYUNIT,*) 'FINAL ITYPE (', ITYPE, ') LIST:', &759:   !WRITE(MYUNIT,*) 'FINAL ITYPE (', ITYPE, ') LIST:', &
763:   !     (ATOMLISTS(ITYPE,IGROUP,K), K=1,ATOMLISTS(ITYPE,IGROUP,0))760:   !     (ATOMLISTS(ITYPE,IGROUP,K), K=1,ATOMLISTS(ITYPE,IGROUP,0))
764:   !WRITE(MYUNIT,*) 'FINLA JTYPE (', JTYPE, ') LIST:', &761:   !WRITE(MYUNIT,*) 'FINLA JTYPE (', JTYPE, ') LIST:', &
765:   !     (ATOMLISTS(JTYPE,JGROUP,K), K=1,ATOMLISTS(JTYPE,JGROUP,0))762:   !     (ATOMLISTS(JTYPE,JGROUP,K), K=1,ATOMLISTS(JTYPE,JGROUP,0))
766:   RETURN763:   RETURN
767:   !764:   !
768: END SUBROUTINE SWAP_LABELS765: END SUBROUTINE SWAP_LABELS
769: !766: !
770: !=============================================================767: !=============================================================
771: !768: !
772: SUBROUTINE FLIP_LABEL(I,NEWTYPE,NP)769: SUBROUTINE FLIP_LABEL(I,NEWTYPE)
773:   !770:   !
774:   ! Change the current label/type of a particle (identified by its771:   ! Change the current label/type of a particle (identified by its
775:   ! permanent index I) to NEWTYPE.772:   ! permanent index I) to NEWTYPE.
776:   !773:   !
777:   ! USE COMMONS, ONLY : ATOMLISTS, INVATOMLISTS, MYUNIT774:   ! USE COMMONS, ONLY : ATOMLISTS, INVATOMLISTS, MYUNIT
778:   USE COMMONS, ONLY: ATOMLISTS, INVATOMLISTS, NSPECIES, ATMASS, &775:   USE COMMONS, ONLY: ATOMLISTS, INVATOMLISTS, NSPECIES, ATMASS, &
779:        READMASST, SPECMASST, SPECMASS, LABELS776:        READMASST, SPECMASST, SPECMASS
780:   !777:   !
781:   IMPLICIT NONE778:   IMPLICIT NONE
782:   !779:   !
783:   INTEGER, INTENT(IN) :: I, NEWTYPE,NP780:   INTEGER, INTENT(IN) :: I, NEWTYPE
784:   INTEGER :: ITYPE,IPLACE,IGROUP,J781:   INTEGER :: ITYPE,IPLACE,IGROUP,J
785:   !782:   !
786:   ITYPE = INVATOMLISTS(I,1)783:   ITYPE = INVATOMLISTS(I,1)
787:   IGROUP = INVATOMLISTS(I,2)784:   IGROUP = INVATOMLISTS(I,2)
788:   IPLACE = INVATOMLISTS(I,3)785:   IPLACE = INVATOMLISTS(I,3)
789:   !786:   !
790:   LABELS(I,NP) = NEWTYPE 
791:   ! 
792:   IF(IPLACE < ATOMLISTS(ITYPE,IGROUP,0)) THEN787:   IF(IPLACE < ATOMLISTS(ITYPE,IGROUP,0)) THEN
793:      ! Fill the hole at IPLACE with the last entry for ITYPE788:      ! Fill the hole at IPLACE with the last entry for ITYPE
794:      J = ATOMLISTS(ITYPE,IGROUP,ATOMLISTS(ITYPE,IGROUP,0))789:      J = ATOMLISTS(ITYPE,IGROUP,ATOMLISTS(ITYPE,IGROUP,0))
795:      ATOMLISTS(ITYPE,IGROUP,IPLACE) = J790:      ATOMLISTS(ITYPE,IGROUP,IPLACE) = J
796:      INVATOMLISTS(J,3) = IPLACE       791:      INVATOMLISTS(J,3) = IPLACE       
797:   ENDIF792:   ENDIF
798:   ATOMLISTS(ITYPE,IGROUP,ATOMLISTS(ITYPE,IGROUP,0)) = 0793:   ATOMLISTS(ITYPE,IGROUP,ATOMLISTS(ITYPE,IGROUP,0)) = 0
799:   ATOMLISTS(ITYPE,IGROUP,0) = ATOMLISTS(ITYPE,IGROUP,0) - 1794:   ATOMLISTS(ITYPE,IGROUP,0) = ATOMLISTS(ITYPE,IGROUP,0) - 1
800:   ! Is there a problem when IPLACE is the last entry, i.e. J = I?795:   ! Is there a problem when IPLACE is the last entry, i.e. J = I?
801:   !796:   !
892:      DO TB=TA+1,NSPECIES(0)887:      DO TB=TA+1,NSPECIES(0)
893:         !888:         !
894:         DO GA = 1,3 ! Double loop over all groups889:         DO GA = 1,3 ! Double loop over all groups
895:            DO GB = 1,3 890:            DO GB = 1,3 
896:               !891:               !
897:               DO IA=1,ATOMLISTS(TA,GA,0) ! Double loop over atoms892:               DO IA=1,ATOMLISTS(TA,GA,0) ! Double loop over atoms
898:                  I=ATOMLISTS(TA,GA,IA)893:                  I=ATOMLISTS(TA,GA,IA)
899:                  DO IB=1,ATOMLISTS(TB,GB,0)894:                  DO IB=1,ATOMLISTS(TB,GB,0)
900:                     J=ATOMLISTS(TB,GB,IB)895:                     J=ATOMLISTS(TB,GB,IB)
901:                     !896:                     !
902:                     CALL SWAP_LABELS(I,J, NP)897:                     CALL SWAP_LABELS(I,J)
903:                     !898:                     !
904:                     NQTOT = NQTOT + 1899:                     NQTOT = NQTOT + 1
905:                     NQ(NP) = NQ(NP) + 1900:                     NQ(NP) = NQ(NP) + 1
906:                     CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN, &901:                     CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN, &
907:                          QDONE,SCREENC)902:                          QDONE,SCREENC)
908:                     IF(QALCSV) CALL PRINT_QUENCH(NP, ITER, '  ')903:                     IF(QALCSV) CALL PRINT_QUENCH(NP, ITER, '  ')
909:                     !904:                     !
910:                     !NTOT = NTOT+1905:                     !NTOT = NTOT+1
911:                     IF(POTEL < E0LO) THEN906:                     IF(POTEL < E0LO) THEN
912:                        NNEG=NNEG+1                       907:                        NNEG=NNEG+1                       
916:                           IF(POTEL < EMIN(1) - ECONV) THEN911:                           IF(POTEL < EMIN(1) - ECONV) THEN
917:                              EMIN(2) = EMIN(1)912:                              EMIN(2) = EMIN(1)
918:                              BESTSWAP(2,1:2) = BESTSWAP(1,1:2)913:                              BESTSWAP(2,1:2) = BESTSWAP(1,1:2)
919:                              EMIN(1) = POTEL914:                              EMIN(1) = POTEL
920:                              BESTSWAP(1,1:2) = (/I,J/)915:                              BESTSWAP(1,1:2) = (/I,J/)
921:                              IF(STEP) XMIN(:) = COORDS(:,NP)916:                              IF(STEP) XMIN(:) = COORDS(:,NP)
922:                           ENDIF917:                           ENDIF
923:                        ENDIF918:                        ENDIF
924:                     ENDIF919:                     ENDIF
925:                     !920:                     !
926:                     CALL SWAP_LABELS(I,J,NP)921:                     CALL SWAP_LABELS(I,J)
927:                     POTEL = E0922:                     POTEL = E0
928:                     COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)923:                     COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)
929:                     !924:                     !
930:                  ENDDO925:                  ENDDO
931:               ENDDO926:               ENDDO
932:               !927:               !
933:            ENDDO928:            ENDDO
934:         ENDDO929:         ENDDO
935:      ENDDO930:      ENDDO
936:   ENDDO931:   ENDDO
952:      ENDIF947:      ENDIF
953:   ELSE948:   ELSE
954:      WRITE(MYUNIT,*)949:      WRITE(MYUNIT,*)
955:   ENDIF950:   ENDIF
956:   !951:   !
957:   IF(STEP) THEN952:   IF(STEP) THEN
958:      IF(NNEG > 0) THEN953:      IF(NNEG > 0) THEN
959:         SEQLENGTH = SEQLENGTH + 1954:         SEQLENGTH = SEQLENGTH + 1
960:         POTEL = EMIN(1)955:         POTEL = EMIN(1)
961:         COORDS(1:3*NATOMS, NP) = XMIN(1:3*NATOMS)956:         COORDS(1:3*NATOMS, NP) = XMIN(1:3*NATOMS)
962:         CALL SWAP_LABELS(BESTSWAP(1,1),BESTSWAP(1,2),NP)957:         CALL SWAP_LABELS(BESTSWAP(1,1),BESTSWAP(1,2))
963:         WRITE(MYUNIT, '(A, G20.10)') &958:         WRITE(MYUNIT, '(A, G20.10)') &
964:              'span_swaps> Executed best swap with dE=', &959:              'span_swaps> Executed best swap with dE=', &
965:              EMIN(1)-E0960:              EMIN(1)-E0
966:      ELSE961:      ELSE
967:         STEP=.FALSE.962:         STEP=.FALSE.
968:      ENDIF963:      ENDIF
969:   ELSEIF(NNEG > 0) THEN964:   ELSEIF(NNEG > 0) THEN
970:      STEP=.TRUE.965:      STEP=.TRUE.
971:   ENDIF966:   ENDIF
972:   !967:   !
1018:   NNEG=0 ! counter for -ve flip gains1013:   NNEG=0 ! counter for -ve flip gains
1019:   !NTOT=01014:   !NTOT=0
1020:   !1015:   !
1021:   DO TA=1,NSPECIES(0) ! loop over species1016:   DO TA=1,NSPECIES(0) ! loop over species
1022:      DO GA = 1,2 ! loop over groups1017:      DO GA = 1,2 ! loop over groups
1023:         DO IA=1,ATOMLISTS(TA,GA,0) ! loop over atoms1018:         DO IA=1,ATOMLISTS(TA,GA,0) ! loop over atoms
1024:            I=ATOMLISTS(TA,GA,IA) ! actual atom index1019:            I=ATOMLISTS(TA,GA,IA) ! actual atom index
1025:            DO TB=1,NSPECIES(0)1020:            DO TB=1,NSPECIES(0)
1026:               IF(TB.NE.TA) THEN ! attempt flip1021:               IF(TB.NE.TA) THEN ! attempt flip
1027:                  !1022:                  !
1028:                  CALL FLIP_LABEL(I,TB, NP)1023:                  CALL FLIP_LABEL(I,TB)
1029:                  !1024:                  !
1030:                  NQTOT = NQTOT + 11025:                  NQTOT = NQTOT + 1
1031:                  NQ(NP) = NQ(NP) + 11026:                  NQ(NP) = NQ(NP) + 1
1032:                  CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN, &1027:                  CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN, &
1033:                       QDONE,SCREENC)1028:                       QDONE,SCREENC)
1034:                  IF(QALCSV) CALL PRINT_QUENCH(NP, ITER, '  ')1029:                  IF(QALCSV) CALL PRINT_QUENCH(NP, ITER, '  ')
1035:                  !1030:                  !
1036:                  !NTOT = NTOT+11031:                  !NTOT = NTOT+1
1037:                  IF(POTEL < E0LO) THEN1032:                  IF(POTEL < E0LO) THEN
1038:                     NNEG=NNEG+1                       1033:                     NNEG=NNEG+1                       
1042:                        IF(POTEL < EMIN(1) - ECONV) THEN1037:                        IF(POTEL < EMIN(1) - ECONV) THEN
1043:                           EMIN(2) = EMIN(1)1038:                           EMIN(2) = EMIN(1)
1044:                           BESTFLIP(2,1:2) = BESTFLIP(1,1:2)1039:                           BESTFLIP(2,1:2) = BESTFLIP(1,1:2)
1045:                           EMIN(1) = POTEL1040:                           EMIN(1) = POTEL
1046:                           BESTFLIP(1,1:2) = (/I,TB/)1041:                           BESTFLIP(1,1:2) = (/I,TB/)
1047:                           IF(STEP) XMIN(:) = COORDS(:,NP)1042:                           IF(STEP) XMIN(:) = COORDS(:,NP)
1048:                        ENDIF1043:                        ENDIF
1049:                     ENDIF1044:                     ENDIF
1050:                  ENDIF1045:                  ENDIF
1051:                  !1046:                  !
1052:                  CALL FLIP_LABEL(I,TA, NP)1047:                  CALL FLIP_LABEL(I,TA)
1053:                  POTEL = E01048:                  POTEL = E0
1054:                  COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)1049:                  COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)
1055:                  !1050:                  !
1056:               ENDIF1051:               ENDIF
1057:            ENDDO1052:            ENDDO
1058:         ENDDO1053:         ENDDO
1059:         !1054:         !
1060:      ENDDO1055:      ENDDO
1061:   ENDDO1056:   ENDDO
1062:   !1057:   !
1077:   !   ENDIF1072:   !   ENDIF
1078:   ELSE1073:   ELSE
1079:      WRITE(MYUNIT,*)1074:      WRITE(MYUNIT,*)
1080:   ENDIF1075:   ENDIF
1081:   !1076:   !
1082:   IF(STEP) THEN1077:   IF(STEP) THEN
1083:      IF(NNEG > 0) THEN1078:      IF(NNEG > 0) THEN
1084:         SEQLENGTH = SEQLENGTH + 11079:         SEQLENGTH = SEQLENGTH + 1
1085:         POTEL = EMIN(1)1080:         POTEL = EMIN(1)
1086:         COORDS(1:3*NATOMS, NP) = XMIN(1:3*NATOMS)1081:         COORDS(1:3*NATOMS, NP) = XMIN(1:3*NATOMS)
1087:         CALL FLIP_LABEL(BESTFLIP(1,1),BESTFLIP(1,2),NP)1082:         CALL FLIP_LABEL(BESTFLIP(1,1),BESTFLIP(1,2))
1088:         WRITE(MYUNIT, '(A, G20.10)') &1083:         WRITE(MYUNIT, '(A, G20.10)') &
1089:              'span_flips> Executed best flip with dE=', &1084:              'span_flips> Executed best flip with dE=', &
1090:              EMIN(1)-E01085:              EMIN(1)-E0
1091:      ELSE1086:      ELSE
1092:         STEP=.FALSE.1087:         STEP=.FALSE.
1093:      ENDIF1088:      ENDIF
1094:   ELSEIF(NNEG > 0) THEN1089:   ELSEIF(NNEG > 0) THEN
1095:      STEP=.TRUE.1090:      STEP=.TRUE.
1096:   ENDIF1091:   ENDIF
1097:   !1092:   !
1220:                  I=ATOMLISTS(TA,GA,IA)1215:                  I=ATOMLISTS(TA,GA,IA)
1221:                  DO IB=1,ATOMLISTS(TB,GB,0)1216:                  DO IB=1,ATOMLISTS(TB,GB,0)
1222:                     J=ATOMLISTS(TB,GB,IB)1217:                     J=ATOMLISTS(TB,GB,IB)
1223:                     !1218:                     !
1224:                     N = N + 11219:                     N = N + 1
1225:                     SWAPS(N,1) = I1220:                     SWAPS(N,1) = I
1226:                     SWAPS(N,2) = J1221:                     SWAPS(N,2) = J
1227:                     !1222:                     !
1228:                     IF(QALCSMODE == 5) THEN1223:                     IF(QALCSMODE == 5) THEN
1229:                        !1224:                        !
1230:                        CALL SWAP_LABELS(I,J,NP)1225:                        CALL SWAP_LABELS(I,J)
1231:                        CALL POTENTIAL(COORDS(:,NP),X0(:),EST(N),.FALSE.,.FALSE.)1226:                        CALL POTENTIAL(COORDS(:,NP),X0(:),EST(N),.FALSE.,.FALSE.)
1232:                        CALL SWAP_LABELS(I,J,NP)1227:                        CALL SWAP_LABELS(I,J)
1233:                        !1228:                        !
1234:                        SORT_EST2: DO K=N,2,-11229:                        SORT_EST2: DO K=N,2,-1
1235:                           IF(EST(K) < EST(K-1)) THEN1230:                           IF(EST(K) < EST(K-1)) THEN
1236:                              E0=EST(K); EST(K)=EST(K-1); EST(K-1)=E01231:                              E0=EST(K); EST(K)=EST(K-1); EST(K-1)=E0
1237:                              I_AB(1:2) = SWAPS(K,1:2)1232:                              I_AB(1:2) = SWAPS(K,1:2)
1238:                              SWAPS(K,1:2) = SWAPS(K-1,1:2)1233:                              SWAPS(K,1:2) = SWAPS(K-1,1:2)
1239:                              SWAPS(K-1,1:2) = I_AB(1:2)1234:                              SWAPS(K-1,1:2) = I_AB(1:2)
1240:                           ELSE1235:                           ELSE
1241:                              EXIT SORT_EST21236:                              EXIT SORT_EST2
1242:                           ENDIF1237:                           ENDIF
1270:   ENDIF1265:   ENDIF
1271:   !1266:   !
1272:   DO I=1,NTOT1267:   DO I=1,NTOT
1273:      !1268:      !
1274:      IF(QALCSMODE==4) THEN    1269:      IF(QALCSMODE==4) THEN    
1275:         J = INT(DPRAND()*DBLE(N)) + 11270:         J = INT(DPRAND()*DBLE(N)) + 1
1276:      ELSEIF(QALCSMODE==5) THEN1271:      ELSEIF(QALCSMODE==5) THEN
1277:         J=I1272:         J=I
1278:      ENDIF1273:      ENDIF
1279:      !1274:      !
1280:      CALL SWAP_LABELS(SWAPS(J,1),SWAPS(J,2),NP)1275:      CALL SWAP_LABELS(SWAPS(J,1),SWAPS(J,2))
1281:      !                                                            1276:      !                                                            
1282:      NQTOT = NQTOT + 11277:      NQTOT = NQTOT + 1
1283:      NQ(NP) = NQ(NP) + 11278:      NQ(NP) = NQ(NP) + 1
1284:      CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)1279:      CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)
1285:      IF(QALCSV) THEN1280:      IF(QALCSV) THEN
1286:         CALL PRINT_QUENCH(NP, ITER, '  ')1281:         CALL PRINT_QUENCH(NP, ITER, '  ')
1287:         IF(QALCSMODE==5) THEN1282:         IF(QALCSMODE==5) THEN
1288:            WRITE(MYUNIT,'(2(A,F15.10))') &1283:            WRITE(MYUNIT,'(2(A,F15.10))') &
1289:                 'scan_nbrhd> dE*= ',EST(J)-E0,' dE= ',POTEL-E0 1284:                 'scan_nbrhd> dE*= ',EST(J)-E0,' dE= ',POTEL-E0 
1290:         ENDIF1285:         ENDIF
1293:      IF(POTEL < E0 - ECONV) THEN1288:      IF(POTEL < E0 - ECONV) THEN
1294:         SEQLENGTH = SEQLENGTH + 11289:         SEQLENGTH = SEQLENGTH + 1
1295:         WRITE(MYUNIT,'(A,I4,A,I4,A,I4,A,F15.10)') &1290:         WRITE(MYUNIT,'(A,I4,A,I4,A,I4,A,F15.10)') &
1296:              'scan_nbrhd> Try ',I,' : ',&1291:              'scan_nbrhd> Try ',I,' : ',&
1297:              SWAPS(J,1),' <-> ',SWAPS(J,2),&1292:              SWAPS(J,1),' <-> ',SWAPS(J,2),&
1298:              ' => dE= ', POTEL-E01293:              ' => dE= ', POTEL-E0
1299:         STEP = .TRUE.1294:         STEP = .TRUE.
1300:         RETURN1295:         RETURN
1301:      ELSE 1296:      ELSE 
1302:         ! Undo the swap and restore the configuration1297:         ! Undo the swap and restore the configuration
1303:         CALL SWAP_LABELS(SWAPS(J,1),SWAPS(J,2),NP)1298:         CALL SWAP_LABELS(SWAPS(J,1),SWAPS(J,2))
1304:         COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)1299:         COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)
1305:         POTEL = E01300:         POTEL = E0
1306:         !1301:         !
1307:         IF(QALCSMODE==4) THEN1302:         IF(QALCSMODE==4) THEN
1308:            ! Remove the rejected swap from list1303:            ! Remove the rejected swap from list
1309:            SWAPS(J,1:2) = SWAPS(N,1:2)1304:            SWAPS(J,1:2) = SWAPS(N,1:2)
1310:            SWAPS(N,1:2) = 01305:            SWAPS(N,1:2) = 0
1311:            N = N - 11306:            N = N - 1
1312:         ENDIF1307:         ENDIF
1313:         !1308:         !
1363:      DO GA = 1,2 ! loop over groups 1 and 21358:      DO GA = 1,2 ! loop over groups 1 and 2
1364:         DO IA=1,ATOMLISTS(TA,GA,0) ! loop over atoms1359:         DO IA=1,ATOMLISTS(TA,GA,0) ! loop over atoms
1365:            I=ATOMLISTS(TA,GA,IA)1360:            I=ATOMLISTS(TA,GA,IA)
1366:            DO TB=1,NSPECIES(0) ! loop over species1361:            DO TB=1,NSPECIES(0) ! loop over species
1367:               IF(TB.NE.TA) THEN1362:               IF(TB.NE.TA) THEN
1368:                  N = N + 11363:                  N = N + 1
1369:                  FLIPS(N,1) = I; FLIPS(N,2) = TB                    1364:                  FLIPS(N,1) = I; FLIPS(N,2) = TB                    
1370:                  !1365:                  !
1371:                  IF(QALCSMODE == 8) THEN1366:                  IF(QALCSMODE == 8) THEN
1372:                     !1367:                     !
1373:                     CALL FLIP_LABEL(I,TB,NP)1368:                     CALL FLIP_LABEL(I,TB)
1374:                     CALL POTENTIAL(COORDS(:,NP),X0(:),EST(N),.FALSE.,.FALSE.)1369:                     CALL POTENTIAL(COORDS(:,NP),X0(:),EST(N),.FALSE.,.FALSE.)
1375:                     IF (SEMIGRAND_MUT) THEN1370:                     IF (SEMIGRAND_MUT) THEN
1376:                        DUMMY=0.0D01371:                        DUMMY=0.0D0
1377:                        DO J=2, NSPECIES(0)1372:                        DO J=2, NSPECIES(0)
1378:                           DUMMY = DUMMY + NSPECIES(J)*SEMIGRAND_MU(J)1373:                           DUMMY = DUMMY + NSPECIES(J)*SEMIGRAND_MU(J)
1379:                        ENDDO1374:                        ENDDO
1380:                        EST(N) = EST(N) - DUMMY1375:                        EST(N) = EST(N) - DUMMY
1381:                     ENDIF1376:                     ENDIF
1382:                     CALL FLIP_LABEL(I,TA,NP)1377:                     CALL FLIP_LABEL(I,TA)
1383:                     !1378:                     !
1384:                     SORT_EST3: DO K=N,2,-11379:                     SORT_EST3: DO K=N,2,-1
1385:                        IF(EST(K) < EST(K-1)) THEN1380:                        IF(EST(K) < EST(K-1)) THEN
1386:                           E0=EST(K); EST(K)=EST(K-1); EST(K-1)=E01381:                           E0=EST(K); EST(K)=EST(K-1); EST(K-1)=E0
1387:                           I_AB(1:2) = FLIPS(K,1:2)1382:                           I_AB(1:2) = FLIPS(K,1:2)
1388:                           FLIPS(K,1:2) = FLIPS(K-1,1:2)1383:                           FLIPS(K,1:2) = FLIPS(K-1,1:2)
1389:                           FLIPS(K-1,1:2) = I_AB(1:2)1384:                           FLIPS(K-1,1:2) = I_AB(1:2)
1390:                        ELSE1385:                        ELSE
1391:                           EXIT SORT_EST31386:                           EXIT SORT_EST3
1392:                        ENDIF1387:                        ENDIF
1419:   DO I=1,NTOT1414:   DO I=1,NTOT
1420:      !1415:      !
1421:      IF(QALCSMODE==7) THEN    1416:      IF(QALCSMODE==7) THEN    
1422:         J = INT(DPRAND()*DBLE(N)) + 11417:         J = INT(DPRAND()*DBLE(N)) + 1
1423:      ELSEIF(QALCSMODE==8) THEN1418:      ELSEIF(QALCSMODE==8) THEN
1424:         J=I1419:         J=I
1425:      ENDIF1420:      ENDIF
1426:      !1421:      !
1427:      TA=INVATOMLISTS(FLIPS(J,1),1) ! Current label1422:      TA=INVATOMLISTS(FLIPS(J,1),1) ! Current label
1428:      TB=FLIPS(J,2) ! new label1423:      TB=FLIPS(J,2) ! new label
1429:      CALL FLIP_LABEL(FLIPS(J,1),TB,NP)1424:      CALL FLIP_LABEL(FLIPS(J,1),TB)
1430:      !                                                            1425:      !                                                            
1431:      NQTOT = NQTOT + 11426:      NQTOT = NQTOT + 1
1432:      NQ(NP) = NQ(NP) + 11427:      NQ(NP) = NQ(NP) + 1
1433:      CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)1428:      CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC)
1434:      IF(QALCSV) THEN1429:      IF(QALCSV) THEN
1435:         CALL PRINT_QUENCH(NP, ITER, '  ')1430:         CALL PRINT_QUENCH(NP, ITER, '  ')
1436:         IF(QALCSMODE==8) THEN1431:         IF(QALCSMODE==8) THEN
1437:            WRITE(MYUNIT,'(2(A,F15.10))') &1432:            WRITE(MYUNIT,'(2(A,F15.10))') &
1438:                 'scan_flip_nbrhd> dE*= ',EST(J)-E0,' dE= ',POTEL-E0 1433:                 'scan_flip_nbrhd> dE*= ',EST(J)-E0,' dE= ',POTEL-E0 
1439:         ENDIF1434:         ENDIF
1442:      IF(POTEL < E0 - ECONV) THEN1437:      IF(POTEL < E0 - ECONV) THEN
1443:         SEQLENGTH = SEQLENGTH + 11438:         SEQLENGTH = SEQLENGTH + 1
1444:         WRITE(MYUNIT,'(A,I4,A,I4,A,I4,A,F15.10)') &1439:         WRITE(MYUNIT,'(A,I4,A,I4,A,I4,A,F15.10)') &
1445:              'scan_flip_nbrhd> Trial ',I,' : ',&1440:              'scan_flip_nbrhd> Trial ',I,' : ',&
1446:              FLIPS(J,1),' ~> ',FLIPS(J,2),&1441:              FLIPS(J,1),' ~> ',FLIPS(J,2),&
1447:              ' => dE= ', POTEL-E01442:              ' => dE= ', POTEL-E0
1448:         STEP = .TRUE.1443:         STEP = .TRUE.
1449:         RETURN1444:         RETURN
1450:      ELSE 1445:      ELSE 
1451:         ! Undo the flip and restore the configuration1446:         ! Undo the flip and restore the configuration
1452:         CALL FLIP_LABEL(FLIPS(J,1),TA,NP)1447:         CALL FLIP_LABEL(FLIPS(J,1),TA)
1453:         COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)1448:         COORDS(1:3*NATOMS, NP) = X0(1:3*NATOMS)
1454:         POTEL = E01449:         POTEL = E0
1455:         !1450:         !
1456:         IF(QALCSMODE==7) THEN1451:         IF(QALCSMODE==7) THEN
1457:            ! Remove the rejected swap from list1452:            ! Remove the rejected swap from list
1458:            FLIPS(J,1:2) = FLIPS(N,1:2)1453:            FLIPS(J,1:2) = FLIPS(N,1:2)
1459:            FLIPS(N,1:2) = 01454:            FLIPS(N,1:2) = 0
1460:            N = N - 11455:            N = N - 1
1461:         ENDIF1456:         ENDIF
1462:         !1457:         !


r31119/quench.F 2016-09-06 15:30:39.622453956 +0100 r31118/quench.F 2016-09-06 15:30:43.018499284 +0100
784:          IF (GOODSTRUCTURE .AND. PERCT .AND. (.NOT. SAVEMULTIMINONLY) .AND. (.NOT. GRADPROBLEMT)) THEN784:          IF (GOODSTRUCTURE .AND. PERCT .AND. (.NOT. SAVEMULTIMINONLY) .AND. (.NOT. GRADPROBLEMT)) THEN
785: !fh301>{{{785: !fh301>{{{
786: !           IF (CHEMSHIFT2) THEN786: !           IF (CHEMSHIFT2) THEN
787: !             IF (DABS(ENERGYCAMSHIFT).GE.0.0001) THEN787: !             IF (DABS(ENERGYCAMSHIFT).GE.0.0001) THEN
788: !               CALL GSAVEIT(POTEL,P,NP)788: !               CALL GSAVEIT(POTEL,P,NP)
789: !               IF (MONITORT) CALL MSAVEIT(POTEL,P,NP)789: !               IF (MONITORT) CALL MSAVEIT(POTEL,P,NP)
790: !             ENDIF790: !             ENDIF
791: !           ELSE791: !           ELSE
792: !fh301>}}}792: !fh301>}}}
793:              IF(NSPECIES(0) > 1) THEN793:              IF(NSPECIES(0) > 1) THEN
794:                 CALL GSAVEIT_MC(EREAL,P,LABELS(1:NATOMS,NP), NP)794:                 CALL GSAVEIT_MC(EREAL,P,INVATOMLISTS(1:NATOMS,1), NP)
795:              ELSE795:              ELSE
796:                 CALL GSAVEIT(EREAL,P,NP)796:                 CALL GSAVEIT(EREAL,P,NP)
797:              ENDIF797:              ENDIF
798:              798:              
799:              IF (NSAVE==0) QMIN(1)=min(QMIN(1),EREAL)799:              IF (NSAVE==0) QMIN(1)=min(QMIN(1),EREAL)
800:              IF (MONITORT) CALL MSAVEIT(EREAL,P,NP)800:              IF (MONITORT) CALL MSAVEIT(EREAL,P,NP)
801: !fh301>{{{801: !fh301>{{{
802: !           ENDIF802: !           ENDIF
803: !fh301>}}}803: !fh301>}}}
804:          ENDIF804:          ENDIF


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0