hdiff output

r31417/mc_gbh.F90 2016-11-01 15:30:10.707736698 +0000 r31416/mc_gbh.F90 2016-11-01 15:30:11.751750791 +0000
131:            !131:            !
132:            ! Accumulate POTEL_LIST with randomly perturbed 132:            ! Accumulate POTEL_LIST with randomly perturbed 
133:            ! energies and check for duplicates133:            ! energies and check for duplicates
134:            POTEL_LIST(IPROC+1) = POTEL134:            POTEL_LIST(IPROC+1) = POTEL
135:            DUPLICATE=.FALSE.135:            DUPLICATE=.FALSE.
136:            IF(IPROC>0) THEN136:            IF(IPROC>0) THEN
137:               DUPLICATE=.FALSE.137:               DUPLICATE=.FALSE.
138:               check_dupes: DO I=IPROC,1,-1138:               check_dupes: DO I=IPROC,1,-1
139:                  IF(ABS(POTEL - POTEL_LIST(I)) < ECONV) THEN139:                  IF(ABS(POTEL - POTEL_LIST(I)) < ECONV) THEN
140:                     DUPLICATE=.TRUE.140:                     DUPLICATE=.TRUE.
 141:                     WRITE(MYUNIT,'(A)') &
 142:                          'mc_gbh> Duplicate skipped.'
141:                     EXIT check_dupes143:                     EXIT check_dupes
142:                  ENDIF144:                  ENDIF
143:               ENDDO check_dupes145:               ENDDO check_dupes
144:            ENDIF146:            ENDIF
145:            !147:            !
146:            IF(.NOT.DUPLICATE .AND. &148:            IF(.NOT.DUPLICATE .AND. &
147:                 ABS(POTEL - E0) > ECONV) THEN149:                 ABS(POTEL - E0) > ECONV) THEN
148:               ! Do surface refinement (parallelised!)150:               ! Do surface refinement (parallelised!)
149:               CALL QALCS_SURF(1, ITERNS, TIME, BRUN, QDONE, SCREENC)151:               CALL QALCS_SURF(1, ITERNS, TIME, BRUN, QDONE, SCREENC)
150:               ! Update back-up on IPROC152:               ! Update back-up on IPROC
151:               IF(MYNODE == IPROC) THEN153:               IF(MYNODE == IPROC) THEN
152:                  XDUM(1:3*NATOMSALLOC) = COORDS(1:3*NATOMSALLOC,1)154:                  XDUM(1:3*NATOMSALLOC) = COORDS(1:3*NATOMSALLOC,1)
153:                  LDUM(1:NATOMSALLOC) = LABELS(1:NATOMSALLOC,1)155:                  LDUM(1:NATOMSALLOC) = LABELS(1:NATOMSALLOC,1)
154:                  EDUM = POTEL156:                  EDUM = POTEL
155:               ENDIF157:               ENDIF
156:            ELSE ! Just inflate energy to prohibit selection158:            ELSE ! Just inflate energy to prohibit selection
157:               WRITE(MYUNIT,'(A)') & 
158:                    'mc_gbh> Duplicate skipped.' 
159:               IF(MYNODE == IPROC) EDUM = 1.0D+99 159:               IF(MYNODE == IPROC) EDUM = 1.0D+99 
160:            ENDIF160:            ENDIF
161:            !161:            !
162:         ENDDO proc_loop162:         ENDDO proc_loop
163:         !163:         !
164:         ! Reset to back-up on each PROC164:         ! Reset to back-up on each PROC
165:         COORDS(1:3*NATOMSALLOC,1) = XDUM(1:3*NATOMSALLOC)165:         COORDS(1:3*NATOMSALLOC,1) = XDUM(1:3*NATOMSALLOC)
166:         LABELS(1:NATOMSALLOC,1) = LDUM(1:NATOMSALLOC)166:         LABELS(1:NATOMSALLOC,1) = LDUM(1:NATOMSALLOC)
167:         POTEL = EDUM167:         POTEL = EDUM
168:         !168:         !
182:         IF (TARGET) IPROCLO=MINLOC(POTEL_LIST(1:NPAR_GBH), 1)-1182:         IF (TARGET) IPROCLO=MINLOC(POTEL_LIST(1:NPAR_GBH), 1)-1
183:         !183:         !
184:         ! Choose an IPROC.184:         ! Choose an IPROC.
185:         CALL CHOOSE_FROM_LIST(NPAR_GBH+1,POTEL_LIST,IPROC)185:         CALL CHOOSE_FROM_LIST(NPAR_GBH+1,POTEL_LIST,IPROC)
186:         !186:         !
187:         IF(IPROC>NPAR_GBH) THEN ! Reject all187:         IF(IPROC>NPAR_GBH) THEN ! Reject all
188:            COORDS(1:3*NATOMSALLOC,1) = X0(1:3*NATOMSALLOC)188:            COORDS(1:3*NATOMSALLOC,1) = X0(1:3*NATOMSALLOC)
189:            LABELS(1:NATOMSALLOC,1) = L0(1:NATOMSALLOC) 189:            LABELS(1:NATOMSALLOC,1) = L0(1:NATOMSALLOC) 
190:            POTEL = E0190:            POTEL = E0
191:            IPROC = 0191:            IPROC = 0
192:            WRITE(MYUNIT, '(A)') 'mc_gbh> Rejected trial(s)'192:            WRITE(MYUNIT, '(A)') 'mc_gbh> Rejected all trial states!'
193:         ELSE193:         ELSE
194:            WRITE(MYUNIT, '(A,I3)') &194:            WRITE(MYUNIT, '(A,I3)') &
195:                 'mc_gbh> Accepted IPROC= ', IPROC    195:                 'mc_gbh> Accepted IPROC= ', IPROC    
196:            IPROC=IPROC-1196:            IPROC=IPROC-1
197:         ENDIF197:         ENDIF
198:         !198:         !
199:      ENDIF199:      ENDIF
200:      !200:      !
201: #ifdef MPI 201: #ifdef MPI 
202:      ! Broadcast IPROC from master202:      ! Broadcast IPROC from master
214:      ! Broadcast HIT from IPROCLO if targetting...214:      ! Broadcast HIT from IPROCLO if targetting...
215:      IF(TARGET) THEN215:      IF(TARGET) THEN
216:         ! Broadcast IPROCLO from master216:         ! Broadcast IPROCLO from master
217:         CALL MPI_BCAST(IPROCLO,1,MPI_INTEGER,&217:         CALL MPI_BCAST(IPROCLO,1,MPI_INTEGER,&
218:              0,MPI_COMM_WORLD,MPIERR)218:              0,MPI_COMM_WORLD,MPIERR)
219:         CALL MPI_BCAST(HIT,1,MPI_LOGICAL,&219:         CALL MPI_BCAST(HIT,1,MPI_LOGICAL,&
220:              IPROCLO,MPI_COMM_WORLD,MPIERR)        220:              IPROCLO,MPI_COMM_WORLD,MPIERR)        
221:      ENDIF221:      ENDIF
222: #endif                     222: #endif                     
223:      !223:      !
224:      !IF(POTEL < E0 - ECONV) THEN224:      IF(POTEL < E0 - ECONV) THEN
225:      IF(ABS(POTEL-E0) > ECONV) THEN 
226:         NREJSTREAK=0225:         NREJSTREAK=0
227:      ELSE226:      ELSE
228:         NREJSTREAK = NREJSTREAK+1227:         NREJSTREAK = NREJSTREAK+1
229:      ENDIF228:      ENDIF
230:      !229:      !
231:      IF(HIT) THEN230:      IF(HIT) THEN
232:         WRITE(MYUNIT,'(A,I3,A,I6)') &231:         WRITE(MYUNIT,'(A,I3,A,I6)') &
233:              'mc_gbh> Target hit stochastically by node ',IPROCLO+1,&232:              'mc_gbh> Target hit stochastically by node ',IPROCLO+1,&
234:              ' on trial ',ITRIAL233:              ' on trial ',ITRIAL
235:         EXIT gbh_loop234:         EXIT gbh_loop
261:   IF(ALLOCATED(AVOIDLIST)) DEALLOCATE(AVOIDLIST)260:   IF(ALLOCATED(AVOIDLIST)) DEALLOCATE(AVOIDLIST)
262:   !261:   !
263:   RETURN262:   RETURN
264:   !263:   !
265: END SUBROUTINE MC_GBH264: END SUBROUTINE MC_GBH
266: 265: 
267: SUBROUTINE CHOOSE_FROM_LIST(N,VALUES,I)266: SUBROUTINE CHOOSE_FROM_LIST(N,VALUES,I)
268:   !267:   !
269:   ! Choose an element of VALUES with Boltzmann probability268:   ! Choose an element of VALUES with Boltzmann probability
270:   !269:   !
271:   USE COMMONS, ONLY : TEMP, QALCSV, MYUNIT, ECONV270:   USE COMMONS, ONLY : TEMP, QALCSV, MYUNIT
272:   !271:   !
273:   IMPLICIT NONE272:   IMPLICIT NONE
274:   !273:   !
275:   INTEGER, INTENT(IN) :: N274:   INTEGER, INTENT(IN) :: N
276:   DOUBLE PRECISION, INTENT(IN) :: VALUES(N)275:   DOUBLE PRECISION, INTENT(IN) :: VALUES(N)
277:   INTEGER, INTENT(OUT) :: I276:   INTEGER, INTENT(OUT) :: I
278:   !277:   !
279:   LOGICAL :: NOREJECT 
280:   DOUBLE PRECISION :: PSUM(N), X, Y, DPRAND, ELOWEST278:   DOUBLE PRECISION :: PSUM(N), X, Y, DPRAND, ELOWEST
281:   !279:   !
282:   !X=0.0D0 ! initialise total sum280:   !X=0.0D0 ! initialise total sum
283:   !DO I=1,N-1281:   !DO I=1,N-1
284:   !   Y = MIN(DEXP(-(VALUES(I)-VALUES(N))/TEMP(1)), 1.0D0)282:   !   Y = MIN(DEXP(-(VALUES(I)-VALUES(N))/TEMP(1)), 1.0D0)
285:   !   X = X + Y283:   !   X = X + Y
286:   !   PSUM(I) = X ! store partial sum284:   !   PSUM(I) = X ! store partial sum
287:   !ENDDO285:   !ENDDO
288:   !PSUM(N) = DBLE(N-1)286:   !PSUM(N) = DBLE(N-1)
289:   !287:   !
290:   ELOWEST=MINVAL(VALUES)288:   ELOWEST=MINVAL(VALUES)
291:   IF(ABS(ELOWEST - VALUES(N)) < ECONV) THEN 
292:      NOREJECT = .FALSE. 
293:   ELSE 
294:      NOREJECT = .TRUE. 
295:   ENDIF 
296:   ! 
297:   X=0.0D0 ! initialise total sum289:   X=0.0D0 ! initialise total sum
298:   DO I=1,N-1290:   DO I=1,N-1
299:      IF( (NOREJECT .AND. VALUES(I) > VALUES(N)+ECONV) .OR. &291:      IF(VALUES(I) > VALUES(N) .AND. VALUES(N) > ELOWEST) THEN
300:           ABS(VALUES(I) - VALUES(N)) < ECONV) THEN 
301:         Y=0.0D0292:         Y=0.0D0
302:      ELSE293:      ELSE
303:         Y=DEXP(-(VALUES(I)-ELOWEST)/TEMP(1))294:         Y=DEXP(-(VALUES(I)-ELOWEST)/TEMP(1))
304:      ENDIF295:      ENDIF
305:      X = X + Y296:      X = X + Y
306:      PSUM(I) = X ! store partial sum297:      PSUM(I) = X ! store partial sum
307:   ENDDO298:   ENDDO
308:   IF(VALUES(N) > ELOWEST + ECONV) THEN299:   IF(VALUES(N) > ELOWEST) THEN
309:      PSUM(N) = PSUM(N-1) ! Disallow rejection300:      PSUM(N) = PSUM(N-1) ! Disallow rejection
310:   ELSE301:   ELSE
311:      PSUM(N) = DBLE(N-1) ! Allow rejection302:      PSUM(N) = DBLE(N-1) ! Allow rejection
312:   ENDIF303:   ENDIF
313:   !304:   !
314:   X=PSUM(N)*DPRAND()305:   X=PSUM(N)*DPRAND()
315:   !306:   !
316:   I=1307:   I=1
317:   DO WHILE (X > PSUM(I))308:   DO WHILE (X > PSUM(I))
318:      I=I+1309:      I=I+1
323:   !   WRITE(MYUNIT, *) 'choose_from_list> psum=', PSUM314:   !   WRITE(MYUNIT, *) 'choose_from_list> psum=', PSUM
324:   !   WRITE(MYUNIT, *) 'choose_from_list> choice=', I315:   !   WRITE(MYUNIT, *) 'choose_from_list> choice=', I
325:   !ENDIF316:   !ENDIF
326:   !317:   !
327:   RETURN318:   RETURN
328:   !319:   !
329: END SUBROUTINE CHOOSE_FROM_LIST320: END SUBROUTINE CHOOSE_FROM_LIST
330: 321: 
331: SUBROUTINE RANDOM_MOVE(RESET,DIM,N,X)322: SUBROUTINE RANDOM_MOVE(RESET,DIM,N,X)
332:   !323:   !
333:   USE COMMONS, ONLY : STEP, ASTEP, MYUNIT, MIEFT324:   USE COMMONS, ONLY :STEP, ASTEP, MYUNIT
334:   !325:   !
335:   IMPLICIT NONE326:   IMPLICIT NONE
336:   !327:   !
337:   LOGICAL, INTENT(IN) :: RESET328:   LOGICAL, INTENT(IN) :: RESET
338:   INTEGER, INTENT(IN) :: DIM, N329:   INTEGER, INTENT(IN) :: DIM, N
339:   DOUBLE PRECISION, INTENT(INOUT) :: X(DIM*N)330:   DOUBLE PRECISION, INTENT(INOUT) :: X(DIM*N)
340:   !331:   !
341:   INTEGER :: I, J, K, N1, N2332:   INTEGER :: I, J, K, N1, N2
342:   DOUBLE PRECISION :: CNTR(DIM), RMIN, RMAX, RALL(N), R, &333:   DOUBLE PRECISION :: CNTR(DIM), RMIN, RMAX, RALL(N), R, &
343:        THETA, PHI, DPRAND, TWOPI334:        THETA, PHI, DPRAND, TWOPI
364:   !355:   !
365:   DO I=1,N356:   DO I=1,N
366:      J=DIM*(I-1)357:      J=DIM*(I-1)
367:      THETA = ACOS(2*DPRAND()-1) !THETA = DPRAND()*PI358:      THETA = ACOS(2*DPRAND()-1) !THETA = DPRAND()*PI
368:      PHI = DPRAND()*TWOPI359:      PHI = DPRAND()*TWOPI
369:      IF(RESET) THEN ! Reset to random coordinates360:      IF(RESET) THEN ! Reset to random coordinates
370:         !        361:         !        
371:         R=RMAX*DPRAND()**(1.0D0/3.0D0) !R = DPRAND()*RMAX362:         R=RMAX*DPRAND()**(1.0D0/3.0D0) !R = DPRAND()*RMAX
372:         X(J+1) = CNTR(1) + R*DSIN(THETA)*DCOS(PHI)363:         X(J+1) = CNTR(1) + R*DSIN(THETA)*DCOS(PHI)
373:         X(J+2) = CNTR(2) + R*DSIN(THETA)*DSIN(PHI)364:         X(J+2) = CNTR(2) + R*DSIN(THETA)*DSIN(PHI)
374:         IF(DIM>2) THEN365:         IF(DIM>2) X(J+3) = CNTR(3) + R*DCOS(THETA)
375:            IF(MIEFT) THEN !<ds656 specific to z=0 substrate!!! 
376:               X(J+3) = RMAX + R*DCOS(THETA) 
377:            ELSE 
378:               X(J+3) = CNTR(3) + R*DCOS(THETA) 
379:            ENDIF 
380:         ENDIF 
381:         !366:         !
382:      ELSEIF(RALL(I) > RMIN) THEN ! Surface angular move367:      ELSEIF(RALL(I) > RMIN) THEN ! Surface angular move
383:         !368:         !
384:         X(J+1) = CNTR(1) + RMAX*DSIN(THETA)*DCOS(PHI)369:         X(J+1) = CNTR(1) + RMAX*DSIN(THETA)*DCOS(PHI)
385:         X(J+2) = CNTR(2) + RMAX*DSIN(THETA)*DSIN(PHI)370:         X(J+2) = CNTR(2) + RMAX*DSIN(THETA)*DSIN(PHI)
386:         IF(DIM>2) X(J+3) = CNTR(3) + RMAX*DCOS(THETA)371:         IF(DIM>2) X(J+3) = CNTR(3) + RMAX*DCOS(THETA)
387:         N2=N2+1372:         N2=N2+1
388:         !373:         !
389:      ELSE ! Interior shake move374:      ELSE ! Interior shake move
390:         !375:         !


r31417/MGupta.f90 2016-11-01 15:30:10.199729841 +0000 r31416/MGupta.f90 2016-11-01 15:30:11.231743771 +0000
 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, GBHT 28:        INVATOMLISTS, STRESS, MYNODE, LABELS
 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)
 47:        ATT_DUM3, RHO_DUM, RHO4HESS(NATOMS) 47:        ATT_DUM3, RHO_DUM, RHO4HESS(NATOMS)
 48:   ! 48:   !
 49:   ! Do some initialisations 49:   ! Do some initialisations
 50:   RIJ(1:NATOMS,1:NATOMS)=0.0D0 50:   RIJ(1:NATOMS,1:NATOMS)=0.0D0
 51:   REP_IJ(1:NATOMS,1:NATOMS)=0.0D0 51:   REP_IJ(1:NATOMS,1:NATOMS)=0.0D0
 52:   ATT_IJ(1:NATOMS,1:NATOMS)=0.0D0 52:   ATT_IJ(1:NATOMS,1:NATOMS)=0.0D0
 53:   RHO(1:NATOMS)=0.0D0 53:   RHO(1:NATOMS)=0.0D0
 54:   RHO4HESS(1:NATOMS)=0.0D0 54:   RHO4HESS(1:NATOMS)=0.0D0
 55:   VT(1:NATOMS) = 0.0D0 55:   VT(1:NATOMS) = 0.0D0
 56:   ! 56:   !
 57:   ! Atoms labels from global LABELS(1:NATOMS,IP) 57:   IP=MYNODE+1 ! Atoms labels from global LABELS(1:NATOMS,IP)
 58:   IF(GBHT) THEN 
 59:      IP=1 
 60:   ELSE 
 61:      IP=MYNODE+1 
 62:   ENDIF 
 63:   ! 58:   !
 64:   DO J1=1,NATOMS-1 59:   DO J1=1,NATOMS-1
 65:      J13 = 3*(J1-1) 60:      J13 = 3*(J1-1)
 66:      T1=LABELS(J1,IP) 61:      T1=LABELS(J1,IP)
 67:      DO J2=J1+1,NATOMS 62:      DO J2=J1+1,NATOMS
 68:         J23 = 3*(J2-1) 63:         J23 = 3*(J2-1)
 69:         T2=LABELS(J2,IP) 64:         T2=LABELS(J2,IP)
 70:         ! 65:         !
 71:         DIST=0.0D0 66:         DIST=0.0D0
 72:         DO I=1,3 67:         DO I=1,3


r31417/potential.f90 2016-11-01 15:30:10.971740262 +0000 r31416/potential.f90 2016-11-01 15:30:12.051754840 +0000
112: !      WRITE(88,'(3G20.10)') ((QMINP(J2, J1), J1=1, 3*NATOMS), J2=1, NSAVE)112: !      WRITE(88,'(3G20.10)') ((QMINP(J2, J1), J1=1, 3*NATOMS), J2=1, NSAVE)
113: !      WRITE(88,'(G20.10)') (TEMP(J1), J1=1, NPAR)113: !      WRITE(88,'(G20.10)') (TEMP(J1), J1=1, NPAR)
114: !      WRITE(88,'(G20.10)') (ACCRAT(J1), J1=1, NPAR)114: !      WRITE(88,'(G20.10)') (ACCRAT(J1), J1=1, NPAR)
115: !      WRITE(88,'(G20.10)') (STEP(J1), J1=1, NPAR)115: !      WRITE(88,'(G20.10)') (STEP(J1), J1=1, NPAR)
116: !      WRITE(88,'(G20.10)') (ASTEP(J1), J1=1, NPAR)116: !      WRITE(88,'(G20.10)') (ASTEP(J1), J1=1, NPAR)
117: !      WRITE(88,'(G20.10)') (OSTEP(J1), J1=1, NPAR)117: !      WRITE(88,'(G20.10)') (OSTEP(J1), J1=1, NPAR)
118: !      CALL SYSTEM('rm ssave')118: !      CALL SYSTEM('rm ssave')
119: !      STOP119: !      STOP
120: !   END IF120: !   END IF
121: 121: 
122: ! ds656> For PTMC runs with MIEF and BOXCENTROID, safer to  
123: !        call BOXCENTROID on every potential call...  
124:    IF(BOXCENTROIDT .AND. PTMC) CALL BOXCENTROID(X) 
125: ! 
126: ! Count the number of potential calls.122: ! Count the number of potential calls.
127:    NPCALL=NPCALL+1123:    NPCALL=NPCALL+1
128: 10 CONTINUE124: 10 CONTINUE
129: 125: 
130:    IF (MULTIPOTT) THEN126:    IF (MULTIPOTT) THEN
131:        IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN127:        IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
132:            XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)128:            XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
133:            CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)               129:            CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)               
134:        END IF130:        END IF
135: 131: 
1300:          CALL COMPRESS(X, GRAD, EREAL, GRADT)1296:          CALL COMPRESS(X, GRAD, EREAL, GRADT)
1301:       END IF1297:       END IF
1302:    END IF1298:    END IF
1303: 1299: 
1304: ! Apply twist forces, if appropriate1300: ! Apply twist forces, if appropriate
1305:    IF (TWISTT) THEN1301:    IF (TWISTT) THEN
1306:       CALL TWIST(X(1:3*NATOMS), NATOMS, GRAD(1:3*NATOMS), EREAL, GRADT)1302:       CALL TWIST(X(1:3*NATOMS), NATOMS, GRAD(1:3*NATOMS), EREAL, GRADT)
1307:    END IF1303:    END IF
1308: 1304: 
1309: ! ds656> Apply a substrate field for keyword MIE_FIELD1305: ! ds656> Apply a substrate field for keyword MIE_FIELD
1310:    !write(myunit,*) 'potential> current1 EREAL=', EREAL 
1311:    IF(MIEFT) CALL MIEF(X,GRAD,EREAL,GRADT,SECT,.FALSE.)1306:    IF(MIEFT) CALL MIEF(X,GRAD,EREAL,GRADT,SECT,.FALSE.)
1312:    !write(myunit,*) 'potential> current2 EREAL=', EREAL1307: 
1313:    !stop 
1314: !1308: !
1315: ! js850> apply harmonic field for keyword HARMONICF1309: ! js850> apply harmonic field for keyword HARMONICF
1316: !1310: !
1317:    IF ( HARMONICF ) THEN1311:    IF ( HARMONICF ) THEN
1318:       DO J1=1, NATOMS1312:       DO J1=1, NATOMS
1319:          IF (HARMONICFLIST(J1)) THEN1313:          IF (HARMONICFLIST(J1)) THEN
1320: !            WRITE(*,'(1x, 6F25.16)') X(1), X(2), X(3), X0(1), X0(2), X0(3)1314: !            WRITE(*,'(1x, 6F25.16)') X(1), X(2), X(3), X0(1), X0(2), X0(3)
1321: !            WRITE(*,*) "#V(1) ", GRAD((J1-1)*3+1)1315: !            WRITE(*,*) "#V(1) ", GRAD((J1-1)*3+1)
1322:             CALL HARMONICFIELD( X(3*J1-2), HARMONICR0(3*J1-2), GRAD(3*J1-2), EREAL)1316:             CALL HARMONICFIELD( X(3*J1-2), HARMONICR0(3*J1-2), GRAD(3*J1-2), EREAL)
1323: !            WRITE(*,*) "#V(1) after ", GRAD((J1-1)*3+1)1317: !            WRITE(*,*) "#V(1) after ", GRAD((J1-1)*3+1)


r31417/QALCS_surface.F90 2016-11-01 15:30:10.455733297 +0000 r31416/QALCS_surface.F90 2016-11-01 15:30:11.487747228 +0000
623:         IF(MIEFT .AND. .NOT. BAD) THEN623:         IF(MIEFT .AND. .NOT. BAD) THEN
624:            BADNESS2: DO J2=1,MIEF_NSITES624:            BADNESS2: DO J2=1,MIEF_NSITES
625:               DUMMY = 0.0D0625:               DUMMY = 0.0D0
626:               DO K=1,3626:               DO K=1,3
627:                  DUMMY = DUMMY + &627:                  DUMMY = DUMMY + &
628:                       (MIEF_SITES(J2,K)-VSITE(K))**2628:                       (MIEF_SITES(J2,K)-VSITE(K))**2
629:               ENDDO629:               ENDDO
630:               IF(DUMMY < TOL3) THEN ! overlap630:               IF(DUMMY < TOL3) THEN ! overlap
631:                  BAD = .TRUE.631:                  BAD = .TRUE.
632:                  EXIT BADNESS2632:                  EXIT BADNESS2
633:               !ELSEIF(DUMMY < TOL4) THEN ! add extra weight 633:               ELSEIF(DUMMY < TOL4) THEN ! add extra weight 
634:               !   VSITE_WEIGHT = VSITE_WEIGHT + 1634:                  VSITE_WEIGHT = VSITE_WEIGHT + 1
635:               ENDIF635:               ENDIF
636:            ENDDO BADNESS2636:            ENDDO BADNESS2
637:         ENDIF637:         ENDIF
638:         !638:         !
639:         IF(.NOT. BAD) THEN639:         IF(.NOT. BAD) THEN
640:            IF(NVSITES < NVSITES_MAX) THEN640:            IF(NVSITES < NVSITES_MAX) THEN
641:               NVSITES=NVSITES+1641:               NVSITES=NVSITES+1
642:               VSITES(NVSITES,1:3) = VSITE(1:3)642:               VSITES(NVSITES,1:3) = VSITE(1:3)
643:               !test643:               !test
644:               !WRITE(MYUNIT,'(A,3(1X,F10.5))') 'gen_sites> NEW V-site:',&644:               !WRITE(MYUNIT,'(A,3(1X,F10.5))') 'gen_sites> NEW V-site:',&
699:   ! Sort VSITES by weight from biggest to smallest (CHECK!)699:   ! Sort VSITES by weight from biggest to smallest (CHECK!)
700:   CALL SORT3_V2(NVSITES,NVSITES_MAX,VSITE_WEIGHTS,VSITES)700:   CALL SORT3_V2(NVSITES,NVSITES_MAX,VSITE_WEIGHTS,VSITES)
701:   ! weed out vacancies whose avergage coordination is below 701:   ! weed out vacancies whose avergage coordination is below 
702:   !702:   !
703:   !ds656> testing...703:   !ds656> testing...
704:   !WRITE(MYUNIT,'(A)') 'gen_vsites> Vsites sorted by occupancy:'  704:   !WRITE(MYUNIT,'(A)') 'gen_vsites> Vsites sorted by occupancy:'  
705:   !DO I1=1,NVSITES705:   !DO I1=1,NVSITES
706:   !   WRITE(MYUNIT,'(A,I6,I3,3(1X,F10.5))') 'gen_sites>', &706:   !   WRITE(MYUNIT,'(A,I6,I3,3(1X,F10.5))') 'gen_sites>', &
707:   !        I1, VSITE_WEIGHTS(I1), VSITES(I1,1:3)707:   !        I1, VSITE_WEIGHTS(I1), VSITES(I1,1:3)
708:   !ENDDO708:   !ENDDO
709:   ! write(*,'(I5)') natoms+nvsites709:   write(*,'(I5)') natoms+nvsites
710:   ! write(*,'(A)') 'Atoms and all 1st generation vacancies'710:   write(*,'(A)') 'Atoms and all 1st generation vacancies'
711:   ! DO I0=1,NSPECIES(0) !ds656> Should not exceed 10           711:   DO I0=1,NSPECIES(0) !ds656> Should not exceed 10           
712:   !    DO I1=1,ATOMLISTS(I0,1,0)712:      DO I1=1,ATOMLISTS(I0,1,0)
713:   !       J1=ATOMLISTS(I0,1,I1)713:         J1=ATOMLISTS(I0,1,I1)
714:   !       WRITE(*,'(A,I1,3(1X,F20.10))') 'A',I0,&714:         WRITE(*,'(A,I1,3(1X,F20.10))') 'A',I0,&
715:   !            (COORDS(J2,NP), J2=3*(J1-1)+1,3*J1)715:              (COORDS(J2,NP), J2=3*(J1-1)+1,3*J1)
716:   !    ENDDO716:      ENDDO
717:   ! ENDDO717:   ENDDO
718:   ! DO I1=1,NVSITES718:   DO I1=1,NVSITES
719:   !    WRITE(*,'(I2,3(1X,F20.10))') VSITE_WEIGHTS(I1),&719:      WRITE(*,'(I2,3(1X,F20.10))') VSITE_WEIGHTS(I1),&
720:   !         (VSITES(I1, J2), J2=1,3)720:           (VSITES(I1, J2), J2=1,3)
721:   ! ENDDO721:   ENDDO
722:   ! STOP722:   STOP
723:   !<ds656 ...testing723:   !<ds656 ...testing
724:   !724:   !
725:   I2=NVSITES725:   I2=NVSITES
726:   IF(VSITE_WEIGHTS(1) > 1) THEN ! or 1?726:   IF(VSITE_WEIGHTS(1) > 1) THEN ! or 1?
727:      ! Weed out all but the largest WEIGHTS (highest coordination)727:      ! Weed out all but the largest WEIGHTS (highest coordination)
728:      !I3=MIN(VSITE_WEIGHTS(1),5) ! May be better?!?728:      !I3=MIN(VSITE_WEIGHTS(1),5) ! May be better?!?
729:      I3=VSITE_WEIGHTS(1) 729:      I3=VSITE_WEIGHTS(1) 
730:      WEED: DO I1=I2,1,-1730:      WEED: DO I1=I2,1,-1
731:         IF(VSITE_WEIGHTS(I1) < I3) THEN731:         IF(VSITE_WEIGHTS(I1) < I3) THEN
732:            NVSITES = NVSITES-1732:            NVSITES = NVSITES-1


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0