hdiff output

r31322/mc_gbh.F90 2016-10-16 15:30:08.196818222 +0100 r31321/mc_gbh.F90 2016-10-16 15:30:08.448821639 +0100
 26:   ! 26:   !
 27:   IMPLICIT NONE 27:   IMPLICIT NONE
 28:   ! 28:   !
 29: #ifdef MPI  29: #ifdef MPI 
 30:   INCLUDE 'mpif.h' 30:   INCLUDE 'mpif.h'
 31:   INTEGER MPIERR 31:   INTEGER MPIERR
 32: #endif 32: #endif
 33:   ! 33:   !
 34:   INTEGER, INTENT(IN) :: NSTEPS 34:   INTEGER, INTENT(IN) :: NSTEPS
 35:   ! 35:   !
 36:   LOGICAL :: RESET 36:   LOGICAL :: ACCEPT
 37:   INTEGER :: I, ITRIAL, ITERNS, BRUN, QDONE, NQTOT, IPROC, IPROCLO, & 37:   INTEGER :: I, ITRIAL, ITERNS, BRUN, QDONE, NQTOT, IPROC, IPROCLO, &
 38:        NDUDSTREAK, L0(NATOMSALLOC) 38:        NDUDSTREAK, L0(NATOMSALLOC)
 39:   DOUBLE PRECISION :: TIME, SCREENC(3*NATOMSALLOC), POTEL, R, & 39:   DOUBLE PRECISION :: TIME, SCREENC(3*NATOMSALLOC), POTEL, EXTRA, R, &
 40:        POTEL_LIST(NPAR_GBH+1), DPRAND, EMARK, X0(3*NATOMSALLOC), E0 40:        POTEL_LIST(NPAR_GBH+1), DPRAND, EMARK, X0(3*NATOMSALLOC), E0
 41:   ! 41:   !
 42:   COMMON /MYPOT/ POTEL 42:   COMMON /MYPOT/ POTEL
 43:   ! 43:   !
 44:   WRITE(MYUNIT, '(A)')  'mc_gbh> Calculating initial energy' 44:   WRITE(MYUNIT, '(A)')  'mc_gbh> Calculating initial energy'
 45:   !WRITE(MYUNIT, *)  'mc_gbh> NATOMSALLOC=', NATOMSALLOC 45:   !WRITE(MYUNIT, *)  'mc_gbh> NATOMSALLOC=', NATOMSALLOC
 46:   CALL FLUSH(MYUNIT) 46:   CALL FLUSH(MYUNIT)
 47:   CALL QUENCH(.FALSE.,1,ITERNS,TIME,BRUN,QDONE,SCREENC) 47:   CALL QUENCH(.FALSE.,1,ITERNS,TIME,BRUN,QDONE,SCREENC)
 48:   NQ(1) = 0 48:   NQ(1) = 0
 49:   WRITE(MYUNIT,& 49:   WRITE(MYUNIT,&
 60:   E0 = POTEL 60:   E0 = POTEL
 61:   ! 61:   !
 62:   EMARK=1.0D+99 62:   EMARK=1.0D+99
 63:   NDUDSTREAK=0 63:   NDUDSTREAK=0
 64:   gbh_loop: DO ITRIAL=1,NSTEPS 64:   gbh_loop: DO ITRIAL=1,NSTEPS
 65:      ! 65:      !
 66:      IF(QALCSV) WRITE(MYUNIT, '(A,I10)') & 66:      IF(QALCSV) WRITE(MYUNIT, '(A,I10)') &
 67:           'mc_gbh> Starting iteration ', ITRIAL      67:           'mc_gbh> Starting iteration ', ITRIAL     
 68:      ! 68:      !
 69:      IF(RESTART.AND.NDUDSTREAK >= NRELAX) THEN 69:      IF(RESTART.AND.NDUDSTREAK >= NRELAX) THEN
 70:         RESET=.TRUE. 70:         EXTRA=5.0D0*STEP(1)
 71:         EMARK=1.0D+99 71:         EMARK=1.0D+99
 72:         E0=EMARK ! Guarantee move by inflating current energy 72:         E0=EMARK ! Guarantee move by inflating current energy
 73:         NDUDSTREAK=0 73:         NDUDSTREAK=0
  74:         WRITE(MYUNIT,'(A,I9)') 'mc_gbh> Big shake and E0 reset.'
 74:      ELSE 75:      ELSE
 75:         RESET=.FALSE. 76:         EXTRA=0.0D0
 76:      ENDIF 77:      ENDIF
 77:      ! 78:      !
 78:      ! --- Stochastic cartesian move ------- 79:      ! --- Stochastic cartesian move -------
 79:      !CALL TAKESTEP(1) 80:      !CALL TAKESTEP(1)
 80:      CALL RANDOM_MOVE(RESET,3,NATOMSALLOC,COORDS(:,1))      81:      CALL SHAKE(3*NATOMSALLOC, COORDS(:,1), EXTRA)     
 81:      ! 82:      !
 82:      !IF (RANDMULTIPERMT.AND.MOD(J1,RANDMULTIPERM_STEP)==0) & 83:      !IF (RANDMULTIPERMT.AND.MOD(J1,RANDMULTIPERM_STEP)==0) &
 83:      !     CALL RANDMULTIPERM(1) ! re-write this routine!!!!! 84:      !     CALL RANDMULTIPERM(1) ! re-write this routine!!!!!
 84:      !IF(BOXCENTROIDT) CALL BOXCENTROID(COORDS(:,1)) 85:      !IF(BOXCENTROIDT) CALL BOXCENTROID(COORDS(:,1))
 85:      ! 86:      !
 86:      ! Quench perturbed state 87:      ! Quench perturbed state
 87:      NQ(1) = NQ(1) + 1 88:      NQ(1) = NQ(1) + 1
 88:      CALL QUENCH(.FALSE.,1,ITERNS,TIME,BRUN,QDONE,SCREENC) 89:      CALL QUENCH(.FALSE.,1,ITERNS,TIME,BRUN,QDONE,SCREENC)
 89:      WRITE(MYUNIT,& 90:      WRITE(MYUNIT,&
 90:           '(A,I10,A,G20.10,A,I5,A,G12.5,A,G20.10,A,F10.1)') & 91:           '(A,I10,A,G20.10,A,I5,A,G12.5,A,G20.10,A,F10.1)') &
210:   !IF(QALCSV) THEN211:   !IF(QALCSV) THEN
211:   !   WRITE(MYUNIT, *) 'choose_from_list> vals=', VALUES212:   !   WRITE(MYUNIT, *) 'choose_from_list> vals=', VALUES
212:   !   WRITE(MYUNIT, *) 'choose_from_list> psum=', PSUM213:   !   WRITE(MYUNIT, *) 'choose_from_list> psum=', PSUM
213:   !   WRITE(MYUNIT, *) 'choose_from_list> choice=', I214:   !   WRITE(MYUNIT, *) 'choose_from_list> choice=', I
214:   !ENDIF215:   !ENDIF
215:   !216:   !
216:   RETURN217:   RETURN
217:   !218:   !
218: END SUBROUTINE CHOOSE_FROM_LIST219: END SUBROUTINE CHOOSE_FROM_LIST
219: 220: 
220: SUBROUTINE RANDOM_MOVE(RESET,DIM,N,X)221: SUBROUTINE SHAKE(N,X,EXTRA)
221:   !222:   !
222:   USE COMMONS, ONLY :STEP, ASTEP, MYUNIT223:   USE COMMONS, ONLY : STEP
223:   !224:   !
224:   IMPLICIT NONE225:   IMPLICIT NONE
225:   !226:   !
226:   LOGICAL, INTENT(IN) :: RESET227:   INTEGER, INTENT(IN) :: N
227:   INTEGER, INTENT(IN) :: DIM, N228:   DOUBLE PRECISION, INTENT(INOUT) :: X(N)
228:   DOUBLE PRECISION, INTENT(INOUT) :: X(DIM*N)229:   DOUBLE PRECISION, INTENT(IN) :: EXTRA
229:   ! 
230:   INTEGER :: I, J, K, N1, N2 
231:   DOUBLE PRECISION :: CNTR(DIM), RMIN, RMAX, RALL(N), R, & 
232:        THETA, PHI, DPRAND, TWOPI 
233:   ! 
234:   DOUBLE PRECISION, PARAMETER :: PI=3.141592654D0 
235:   TWOPI=2.0D0*PI 
236:   ! 
237:   CALL FIND_CENTROID(DIM,N,X,CNTR) 
238:   ! 
239:   RMAX=0.0D0 
240:   DO I=1,N 
241:      J=DIM*(I-1) 
242:      R=0.0D0 
243:      DO K=1,DIM 
244:         R = R + (X(I+K)-CNTR(K))**2 
245:      ENDDO 
246:      RALL(I) = DSQRT(R) 
247:      IF(RMAX < RALL(I)) RMAX = RALL(I) 
248:   ENDDO 
249:   ! 
250:   R=(DPRAND()-0.5D0)*2.0D0 ! 50% chance of no angular move 
251:   RMIN = RMAX - R*ASTEP(1) 
252:   N1=0; N2=0 
253:   ! 
254:   DO I=1,N 
255:      J=DIM*(I-1) 
256:      IF(RESET) THEN ! Reset to random coordinates 
257:         ! 
258:         THETA = DPRAND()*PI 
259:         PHI = DPRAND()*TWOPI 
260:         R = DPRAND()*RMAX 
261:         X(J+1) = CNTR(1) + R*DSIN(THETA)*DCOS(PHI) 
262:         X(J+2) = CNTR(2) + R*DSIN(THETA)*DSIN(PHI) 
263:         IF(DIM>2) X(J+3) = CNTR(3) + R*DCOS(THETA) 
264:         ! 
265:      ELSEIF(RALL(I) > RMIN) THEN ! Surface angular move 
266:         ! 
267:         THETA = DPRAND()*PI 
268:         PHI = DPRAND()*TWOPI 
269:         X(J+1) = CNTR(1) + RMAX*DSIN(THETA)*DCOS(PHI) 
270:         X(J+2) = CNTR(2) + RMAX*DSIN(THETA)*DSIN(PHI) 
271:         IF(DIM>2) X(J+3) = CNTR(3) + RMAX*DCOS(THETA) 
272:         N2=N2+1 
273:         ! 
274:      ELSE ! Interior shake move 
275:         ! 
276:         DO K=1,DIM 
277:            R=(DPRAND()-0.5D0)*2.0D0 
278:            X(J+K)=X(J+K) + STEP(1)*R 
279:         ENDDO 
280:         N1=N1+1 
281:         ! 
282:      ENDIF 
283:   ENDDO 
284:   ! 
285:   IF(RESET) THEN 
286:      WRITE(MYUNIT,'(A)') 'mc_gbh> Configuration reset!' 
287:   !ELSEIF(N2>0) THEN 
288:   !   WRITE(MYUNIT, '(A,I6)') 'mc_gbh> Atoms rotated: ', N2 
289:   ENDIF 
290:   ! 
291:   RETURN 
292:   ! 
293: END SUBROUTINE RANDOM_MOVE 
294: ! 
295: SUBROUTINE FIND_CENTROID(DIM,N,X,CNTR) 
296:   ! 
297:   IMPLICIT NONE 
298:   ! 
299:   INTEGER, INTENT(IN) :: DIM, N 
300:   DOUBLE PRECISION, INTENT(IN) :: X(DIM*N) 
301:   DOUBLE PRECISION, INTENT(OUT) :: CNTR(DIM) 
302:   !230:   !
303:   INTEGER :: I, J, K231:   INTEGER :: I
 232:   DOUBLE PRECISION :: DPRAND, R
304:   !233:   !
305:   CNTR(:) = 0.0D0 
306:   DO I=1,N234:   DO I=1,N
307:      J=3*(I-1)235:      R=(DPRAND()-0.5D0)*2.0D0
308:      DO K=1,DIM236:      X(I)=X(I) + STEP(1)*R + EXTRA*SIGN(1.0D0,R)
309:         CNTR(K) = CNTR(K) + X(J+K) 
310:      ENDDO 
311:   ENDDO237:   ENDDO
312:   CNTR(:) = CNTR(:)/DBLE(N) 
313:   !238:   !
314:   RETURN239:   RETURN
315:   !240:   !
316: END SUBROUTINE FIND_CENTROID241: END SUBROUTINE SHAKE


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0