hdiff output

r29933/commons.f90 2016-03-16 18:33:39.879125803 +0000 r29932/commons.f90 2016-03-16 18:33:41.327140692 +0000
 80:      &                 KINT, INTFREEZETOL, IMSEPMIN, IMSEPMAX, CONCUTABS, CONCUTFRAC, & 80:      &                 KINT, INTFREEZETOL, IMSEPMIN, IMSEPMAX, CONCUTABS, CONCUTFRAC, &
 81:      &                 LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, & 81:      &                 LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, &
 82:      &                 INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, & 82:      &                 INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, &
 83:      &                 INTLJTOL, INTLJDEL, INTLJEPS, REPCON, INTDGUESS, CHECKREPCUTOFF, INTMINFAC, FREEZETOL, & 83:      &                 INTLJTOL, INTLJDEL, INTLJEPS, REPCON, INTDGUESS, CHECKREPCUTOFF, INTMINFAC, FREEZETOL, &
 84:      &                 LOCALPERMCUT, LOCALPERMCUT2, INTCONCUT, QCIRADSHIFT, MLPLAMBDA 84:      &                 LOCALPERMCUT, LOCALPERMCUT2, INTCONCUT, QCIRADSHIFT, MLPLAMBDA
 85:  85: 
 86:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, & 86:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, &
 87:      &        OTPT, LJMFT, STRANDT, PAHT, SWT, MSTRANST, STOCKT, STICKYT, BLNT, MYSDT, FREEZERES, CENTXY, & 87:      &        OTPT, LJMFT, STRANDT, PAHT, SWT, MSTRANST, STOCKT, STICKYT, BLNT, MYSDT, FREEZERES, CENTXY, &
 88:      &        MSORIGT, SQUEEZET, PERIODIC, SCT, MSCT, MGUPTAT, RESIZET, TIP, RIGID, CALCQT, MPIT, JMT, LJCOULT, SETCENT, & 88:      &        MSORIGT, SQUEEZET, PERIODIC, SCT, MSCT, MGUPTAT, RESIZET, TIP, RIGID, CALCQT, MPIT, JMT, LJCOULT, SETCENT, &
 89:      &        SORTT, HIT, SAVEQ, PARALLELT, FIXD, RKMIN, BSMIN, PERMDIST, PERMOPT, BSWL, BSPT, BSPTRESTART, & 89:      &        SORTT, HIT, SAVEQ, PARALLELT, FIXD, RKMIN, BSMIN, PERMDIST, PERMOPT, BSWL, BSPT, BSPTRESTART, &
 90:      &        SYMMETRIZE, SYMMETRIZECSM, PRINT_PTGRP, PRINT_MINDATA, DUMPT, NEON, ARGON, P46, NORESET, TABOOT, EVSTEPT, PACHECO, DL_POLY, QUCENTRE, & 90:      &        SYMMETRIZE, SYMMETRIZECSM, PRINT_PTGRP, DUMPT, NEON, ARGON, P46, NORESET, TABOOT, EVSTEPT, PACHECO, DL_POLY, QUCENTRE, &
 91:      &        STAR, PLUS, TWOPLUS, GROUND, DIPOLE, DFTBT, DFTBCT, SW, SUPERSTEP, EAMLJT, PBGLUET, TRACKDATAT, & 91:      &        STAR, PLUS, TWOPLUS, GROUND, DIPOLE, DFTBT, DFTBCT, SW, SUPERSTEP, EAMLJT, PBGLUET, TRACKDATAT, &
 92:      &        EAMALT, ALGLUET, MGGLUET, GUPTAT, LJATT, FST, DECAY, COOP, FIXBIN, GAUSST, QUENCHDOS, FIXDIHEFLAG, & 92:      &        EAMALT, ALGLUET, MGGLUET, GUPTAT, LJATT, FST, DECAY, COOP, FIXBIN, GAUSST, QUENCHDOS, FIXDIHEFLAG, &
 93:      &        FRAUSIT, ANGST, SELFT, STEPOUT, WENZEL, THRESHOLDT, THOMSONT, MULLERBROWNT, CHARMMENERGIES, & 93:      &        FRAUSIT, ANGST, SELFT, STEPOUT, WENZEL, THRESHOLDT, THOMSONT, MULLERBROWNT, CHARMMENERGIES, &
 94:      &        PROJ, RGCL2, TOSI, WELCH, AXTELL, AMBER, FIXIMAGE, BINARY, SHIFTCUT, ARNO, TUNNELT, TWOD, &  94:      &        PROJ, RGCL2, TOSI, WELCH, AXTELL, AMBER, FIXIMAGE, BINARY, SHIFTCUT, ARNO, TUNNELT, TWOD, & 
 95:      &        BLJCLUSTER, BLJCLUSTER_NOCUT, COMPRESST, FIX, FIXT, BFGS, LBFGST, DBRENTT, DZTEST, FNI, FAL, CPMD, TNT, ZETT1, & 95:      &        BLJCLUSTER, BLJCLUSTER_NOCUT, COMPRESST, FIX, FIXT, BFGS, LBFGST, DBRENTT, DZTEST, FNI, FAL, CPMD, TNT, ZETT1, &
 96:      &        ZETT2, RESTART, CONJG, NEWRESTART, AVOID, NATBT, DIFFRACTT, CHRMMT, INTMINT, LB2T, &  96:      &        ZETT2, RESTART, CONJG, NEWRESTART, AVOID, NATBT, DIFFRACTT, CHRMMT, INTMINT, LB2T, & 
 97:      &        PTMC, BINSTRUCTURES, PROGRESS, MODEL1T, NEWRESTART_MD, CHANGE_TEMP, NOCISTRANS, CHECKCHIRALITY, & 97:      &        PTMC, BINSTRUCTURES, PROGRESS, MODEL1T, NEWRESTART_MD, CHANGE_TEMP, NOCISTRANS, CHECKCHIRALITY, &
 98:      &        GBT, GBDT, GBDPT, GEMT, LINRODT, RADIFT, CAPBINT, DBPT, DBPTDT, DMBLMT, DMBLPYT, EFIELDT, PAHAT, STOCKAAT, MORSEDPT, & 98:      &        GBT, GBDT, GBDPT, GEMT, LINRODT, RADIFT, CAPBINT, DBPT, DBPTDT, DMBLMT, DMBLPYT, EFIELDT, PAHAT, STOCKAAT, MORSEDPT, &
 99:      &        MSGBT, MSTBINT, MMRSDPT, MSSTOCKT, LWOTPT, NCAPT, NPAHT, NTIPT, PAHW99T, ELLIPSOIDT, GAYBERNET,&  99:      &        MSGBT, MSTBINT, MMRSDPT, MSSTOCKT, LWOTPT, NCAPT, NPAHT, NTIPT, PAHW99T, ELLIPSOIDT, GAYBERNET,& 
100:      &        MULTPAHAT, PAPT, PAPBINT, PAPJANT, PTSTSTT, SHIFTED, SILANET, TDHDT, DDMT, WATERDCT, WATERKZT, CHECKDT, CHECKMARKOVT,& 100:      &        MULTPAHAT, PAPT, PAPBINT, PAPJANT, PTSTSTT, SHIFTED, SILANET, TDHDT, DDMT, WATERDCT, WATERKZT, CHECKDT, CHECKMARKOVT,& 


r29933/finalio.f90 2016-03-16 18:33:40.075127818 +0000 r29932/finalio.f90 2016-03-16 18:33:41.523142706 +0000
 21:     USE GENRIGID, ONLY : RIGIDINIT 21:     USE GENRIGID, ONLY : RIGIDINIT
 22:     USE MODAMBER 22:     USE MODAMBER
 23:     USE MODAMBER9, ONLY : COORDS1,IH,M04,AMBFINALIO_NODE 23:     USE MODAMBER9, ONLY : COORDS1,IH,M04,AMBFINALIO_NODE
 24:     USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_FINISH, AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, & 24:     USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_FINISH, AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, &
 25:                                       AMBER12_WRITE_XYZ 25:                                       AMBER12_WRITE_XYZ
 26:     USE QMODULE 26:     USE QMODULE
 27:     USE MODCHARMM 27:     USE MODCHARMM
 28:     USE AMHGLOBALS, ONLY:NMRES,IRES 28:     USE AMHGLOBALS, ONLY:NMRES,IRES
 29:     USE BGUPMOD 29:     USE BGUPMOD
 30:     USE PERMU 30:     USE PERMU
 31:     USE MODHESS, ONLY : HESS, MASSWT 
 32:  31: 
 33:     IMPLICIT NONE 32:     IMPLICIT NONE
 34:  33: 
 35:     !   MCP 34:     !   MCP
 36:     INTEGER III, I3,  GLY_COUNT, ID, NUMCRD, NUMPRO, NCPHST, GETUNIT, AMHUNIT1, AMHUNIT2, LUNIT, NSYMOPS 35:     INTEGER III, I3,  GLY_COUNT, ID, NUMCRD, NUMPRO, NCPHST, GETUNIT, AMHUNIT1, AMHUNIT2, LUNIT, NSYMOPS
 37:     INTEGER J1, J2, J3, J4, J5, MYUNIT2, I1, NDUMMY, MYUNIT3, NC, NRBS1, NRBS2, MYUNIT4 36:     INTEGER J1, J2, J3, J4, J5, MYUNIT2, I1, NDUMMY, MYUNIT3, NC, NRBS1, NRBS2, MYUNIT4
 38:     DOUBLE PRECISION RBCOORDS(NRBSITES*3), DCOORDS(3*NATOMS), EDUMMY 37:     DOUBLE PRECISION RBCOORDS(NRBSITES*3), DCOORDS(3*NATOMS), EDUMMY
 39:     DOUBLE PRECISION P3(3,3), P(3), DU(3), RMI(3,3), DRMI(3,3), PI, PHI, THT, CHI 38:     DOUBLE PRECISION P3(3,3), P(3), DU(3), RMI(3,3), DRMI(3,3), PI, PHI, THT, CHI
 40:     !DOUBLE PRECISION, ALLOCATABLE :: XCOORDS(:), YCOORDS(:) 39:     !DOUBLE PRECISION, ALLOCATABLE :: XCOORDS(:), YCOORDS(:)
 41:     CHARACTER(LEN=13) J1CHAR,J1CHAR2                  !for gay-berne output files 40:     CHARACTER(LEN=13) J1CHAR,J1CHAR2                  !for gay-berne output files
 57:     !  AMH 56:     !  AMH
 58:     CHARACTER(LEN=3) :: RES_TYPE 57:     CHARACTER(LEN=3) :: RES_TYPE
 59:     CHARACTER(LEN=2) :: ATOM_TYPE 58:     CHARACTER(LEN=2) :: ATOM_TYPE
 60:     CHARACTER*1 COUNTTT 59:     CHARACTER*1 COUNTTT
 61:     INTEGER COUNTT 60:     INTEGER COUNTT
 62:     DOUBLE PRECISION  PPPCORD(NMRES*3*3,3,3,5) 61:     DOUBLE PRECISION  PPPCORD(NMRES*3*3,3,3,5)
 63:     EXTERNAL NUM_TO_CHAR 62:     EXTERNAL NUM_TO_CHAR
 64:     DOUBLE PRECISION TEND 63:     DOUBLE PRECISION TEND
 65:     INTEGER ITERATIONS, BRUN,QDONE, NP 64:     INTEGER ITERATIONS, BRUN,QDONE, NP
 66:     DOUBLE PRECISION SCREENC(3*NATOMSALLOC), TIME 65:     DOUBLE PRECISION SCREENC(3*NATOMSALLOC), TIME
 67:      
 68:     ! ds656> Extra stuff for HSA in OPTIM's min.data format 
 69:     INTEGER :: NUM_ZERO_EVS, INFO 
 70:     DOUBLE PRECISION :: LOG_PROD, EVALS(3*NATOMS), WORK(9*NATOMS) 
 71:      
 72:  66: 
 73:     PI = 4.D0*DATAN(1.D0) 67:     PI = 4.D0*DATAN(1.D0)
 74:      68:     
 75:     !ds656> test 69:     !ds656> test
 76:     !write(*,*) "finalio> START, PI=", PI 70:     !write(*,*) "finalio> START, PI=", PI
 77:      71:     
 78:     IF (PERMOPT) THEN 72:     IF (PERMOPT) THEN
 79:        NSAVE=2 73:        NSAVE=2
 80:        QMIN(2)=0.0D0 74:        QMIN(2)=0.0D0
 81:        QMINP(2,:)=FIN(:) 75:        QMINP(2,:)=FIN(:)
199:             WRITE(MYUNIT2,*) NATOMS+2193:             WRITE(MYUNIT2,*) NATOMS+2
200:         ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT) THEN194:         ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT) THEN
201:             WRITE(MYUNIT2,*) NATOMS/2195:             WRITE(MYUNIT2,*) NATOMS/2
202:         ELSE IF (AMHT) THEN196:         ELSE IF (AMHT) THEN
203:             WRITE(MYUNIT2,*) NMRES*3197:             WRITE(MYUNIT2,*) NMRES*3
204:         ELSE198:         ELSE
205:             WRITE(MYUNIT2,*) NATOMS199:             WRITE(MYUNIT2,*) NATOMS
206:         ENDIF200:         ENDIF
207:         !        IF (CSMT.AND.DEBUG) WRITE(MYUNIT,'(A,I6,2G20.10)') 'finalio> J1,QMIN,QMINAV=',J1,QMIN(J1),QMINAV(J1)201:         !        IF (CSMT.AND.DEBUG) WRITE(MYUNIT,'(A,I6,2G20.10)') 'finalio> J1,QMIN,QMINAV=',J1,QMIN(J1),QMINAV(J1)
208:         !202:         !
209:         ! ds656> Either append point group or print stuff for HSA203:         ! ds656> print point-group symmetry if required
210:         IF(PRINT_MINDATA.OR.PRINT_PTGRP) THEN204:         IF(PRINT_PTGRP) THEN
211:            ! 
212:            WRITE(MYUNIT,'(A,I4)') 'finalio> Analysing point group of minimum', J1205:            WRITE(MYUNIT,'(A,I4)') 'finalio> Analysing point group of minimum', J1
213:            CM(1:3) = 0.0D0; NSYMOPS=0; SYMOPS(:,:,:) = 0.0D0206:            CM(1:3) = 0.0D0; NSYMOPS=0; SYMOPS(:,:,:) = 0.0D0
214:            !ds656> TEST207:            !ds656> TEST
215:            J3=1 ! Determine majority label/species208:            J3=1 ! Determine majority label/species
216:            DO J2=2,NSPECIES(0)209:            DO J2=2,NSPECIES(0)
217:               IF(NSPECIES(J2) > NSPECIES(J3)) J3 = J2210:               IF(NSPECIES(J2) > NSPECIES(J3)) J3 = J2
218:            ENDDO211:            ENDDO
219:            ! Determine the overall point-group (for ALL species!).212:            ! Determine the overall point-group (for ALL species!).
220:            CALL PGSYM( NATOMS, QMINP(J1,1:3*NATOMS), &213:            CALL PGSYM( NATOMS, QMINP(J1,1:3*NATOMS), &
221:                 QMINT(J1,1:NATOMS), PGSYMTOLS, 2, J3,&214:                 QMINT(J1,1:NATOMS), PGSYMTOLS, 2, J3,&
222:                 CM, NSYMOPS, SYMOPS, FPGRP )215:                 CM, NSYMOPS, SYMOPS, FPGRP )
223:            !216:            WRITE(MYUNIT2,99) J1, QMIN(J1), FF(J1), NPCALL_QMIN(J1), &
224:            IF(PRINT_PTGRP) THEN ! Just append point group217:                 FPGRP, NSYMOPS
225:               WRITE(MYUNIT2,99) J1, QMIN(J1), FF(J1), &218: 99         FORMAT('Energy of minimum ',I6,'=',G20.10, &
226:                    NPCALL_QMIN(J1), FPGRP, NSYMOPS219:                 ' first found at step ',I8,' after ', &
227: 99            FORMAT('Energy of minimum ',I6,'=',G20.10, &220:                 I20,' function calls. Point group ',A4, &
228:                    ' first found at step ',I8,' after ', & 
229:                    I20,' function calls. Point group ',A4, & 
230:                 ' of order ', I3,'.')221:                 ' of order ', I3,'.')
231:            ELSE ! Print stuff for HSA in OPTIM's min.data format 
232:               ! 
233:               NUM_ZERO_EVS=6 
234:               IF (NATOMS.EQ.2) NUM_ZERO_EVS=5 
235:               IF (MIEFT) NUM_ZERO_EVS=0 
236:               ! 
237:               IF (ALLOCATED(HESS)) DEALLOCATE(HESS) 
238:               ALLOCATE(HESS(3*NATOMS,3*NATOMS)) 
239:               IF(NSPECIES(0)>1) & ! ds656> this also sets ATMASS 
240:                    CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1) 
241:               CALL POTENTIAL(QMINP(J1,1:3*NATOMS), EVALS, EDUMMY, & 
242:                    .TRUE., .TRUE.) ! EVALS is dummy for gradient here 
243:               IF(ABS(EDUMMY-QMIN(J1)) > ECONV) THEN 
244:                  WRITE(MYUNIT,'(A)') & 
245:                       'finalio> WARNING: Unexpected energy value!' 
246:               ENDIF 
247:               CALL MASSWT()  
248:               ! 
249:               WRITE(MYUNIT, '(A,I6)') & 
250:                    'finalio> Using DSYEV to calculate Hessian eigenvalues for minimum ', J1 
251:               CALL DSYEV('V','U',3*NATOMS,HESS,3*NATOMS,EVALS,& 
252:                    WORK,9*NATOMS,INFO) 
253:               ! Log prod. of frequencies is log prod. of sqrt of evals  
254:               LOG_PROD = & 
255:                    SUM(DLOG(EVALS(NUM_ZERO_EVS+1:3*NATOMS))) 
256:               WRITE(MYUNIT,'(A,G20.10)') & 
257:                    'finalio> Log product of frequencies= ', & 
258:                    LOG_PROD 
259:               WRITE(MYUNIT2, '(2F20.10,I6)') & 
260:                    EDUMMY, LOG_PROD, NSYMOPS 
261:               ! 
262:            ENDIF 
263:            ! 
264:         ELSE ! <ds656222:         ELSE ! <ds656
265:            WRITE(MYUNIT2,10) J1, QMIN(J1), FF(J1), NPCALL_QMIN(J1)223:            WRITE(MYUNIT2,10) J1, QMIN(J1), FF(J1), NPCALL_QMIN(J1)
266: 10         FORMAT('Energy of minimum ',I6,'=',G20.10, &224: 10         FORMAT('Energy of minimum ',I6,'=',G20.10, &
267:                 ' first found at step ',I8,' after ',I20,' function calls')225:                 ' first found at step ',I8,' after ',I20,' function calls')
268:         ENDIF226:         ENDIF
269:         !227:         !
270:         IF (MSORIGT.OR.FRAUSIT) THEN228:         IF (MSORIGT.OR.FRAUSIT) THEN
271:             WRITE(MYUNIT2,20) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))229:             WRITE(MYUNIT2,20) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))
272: 20          FORMAT('Si',3F20.10)230: 20          FORMAT('Si',3F20.10)
273:         ELSE IF (MSTRANST) THEN231:         ELSE IF (MSTRANST) THEN


r29933/keywords.f 2016-03-16 18:33:40.275129871 +0000 r29932/keywords.f 2016-03-16 18:33:41.719144721 +0000
876: 876: 
877:       UNSTRING='UNDEFINED'877:       UNSTRING='UNDEFINED'
878:       BOXSIZE=20.D0878:       BOXSIZE=20.D0
879:       SPHERERAD=20.D0879:       SPHERERAD=20.D0
880:       NCHENCALLS=0880:       NCHENCALLS=0
881:       NATBT=.FALSE.881:       NATBT=.FALSE.
882:       MAXERISE=1.0D-10882:       MAXERISE=1.0D-10
883:       MAXERISE_SET=.FALSE.883:       MAXERISE_SET=.FALSE.
884:       MAXEFALL=-HUGE(1.0D0)884:       MAXEFALL=-HUGE(1.0D0)
885:       PRINT_PTGRP=.FALSE.885:       PRINT_PTGRP=.FALSE.
886:       PRINT_MINDATA=.FALSE. 
887:       SYMMETRIZE=.FALSE.886:       SYMMETRIZE=.FALSE.
888:       SYMMETRIZECSM=.FALSE.887:       SYMMETRIZECSM=.FALSE.
889:       NSYMINTERVAL=10888:       NSYMINTERVAL=10
890:       SYMTOL1=0.1D0889:       SYMTOL1=0.1D0
891:       SYMTOL2=0.1D0890:       SYMTOL2=0.1D0
892:       SYMTOL3=0.1D0891:       SYMTOL3=0.1D0
893:       SYMTOL4=0.1D0892:       SYMTOL4=0.1D0
894:       SYMTOL5=0.1D0893:       SYMTOL5=0.1D0
895:       PGSYMTOLS(1) = 1.0D-3     ! (Normalised) eigenvalue tolerance894:       PGSYMTOLS(1) = 1.0D-3     ! (Normalised) eigenvalue tolerance
896:       PGSYMTOLS(2) = 2.0D-1     ! Overlap (distance) tolerance895:       PGSYMTOLS(2) = 2.0D-1     ! Overlap (distance) tolerance
5310:          IF(ALLOCATED(SPECMASS)) DEALLOCATE(SPECMASS)5309:          IF(ALLOCATED(SPECMASS)) DEALLOCATE(SPECMASS)
5311:          ALLOCATE(SPECMASS(1:NSPECIES(0)))5310:          ALLOCATE(SPECMASS(1:NSPECIES(0)))
5312:          IF(NITEMS .GT. NSPECIES(0)) THEN5311:          IF(NITEMS .GT. NSPECIES(0)) THEN
5313:             DO J1=1,NSPECIES(0)5312:             DO J1=1,NSPECIES(0)
5314:                CALL READF(SPECMASS(J1))5313:                CALL READF(SPECMASS(J1))
5315:             ENDDO5314:             ENDDO
5316:          ELSE5315:          ELSE
5317:             WRITE(MYUNIT,'(A)')' keywords> ERROR *** not enough arguments for SPECMASS'5316:             WRITE(MYUNIT,'(A)')' keywords> ERROR *** not enough arguments for SPECMASS'
5318:             STOP5317:             STOP
5319:          ENDIF5318:          ENDIF
5320:          IF (ALLOCATED(ATMASS)) DEALLOCATE(ATMASS) 
5321:          ALLOCATE(ATMASS(NATOMSALLOC)) 
5322:          ATMASS(:) = 1.0D0 ! Just initialisation 
5323: !5319: !
5324: ! Low temperature replica will use min.data and points.min entries.5320: ! Low temperature replica will use min.data and points.min entries.
5325: !5321: !
5326:       ELSE IF (WORD.EQ.'RESERVOIR') THEN5322:       ELSE IF (WORD.EQ.'RESERVOIR') THEN
5327:          RESERVOIRT=.TRUE.5323:          RESERVOIRT=.TRUE.
5328:          IF (NITEMS.GT.2) THEN5324:          IF (NITEMS.GT.2) THEN
5329:             CALL READI(USERES)5325:             CALL READI(USERES)
5330:             CALL READF(RES_PSWAP)5326:             CALL READF(RES_PSWAP)
5331:          ENDIF5327:          ENDIF
5332:          INQUIRE(FILE='min.data',EXIST=YESNO)5328:          INQUIRE(FILE='min.data',EXIST=YESNO)
6707:          SETCHIRAL=.TRUE.6703:          SETCHIRAL=.TRUE.
6708:       ELSE IF (WORD.EQ.'SETCHIRALGENERIC') THEN6704:       ELSE IF (WORD.EQ.'SETCHIRALGENERIC') THEN
6709:          SETCHIRALGENERIC=.TRUE.6705:          SETCHIRALGENERIC=.TRUE.
6710: !6706: !
6711:       ELSE IF (WORD .EQ. 'PRINT_PTGRP') THEN6707:       ELSE IF (WORD .EQ. 'PRINT_PTGRP') THEN
6712:          PRINT_PTGRP = .TRUE.6708:          PRINT_PTGRP = .TRUE.
6713:          IF (NITEMS.GT.1) CALL READF(SYMTOL1)6709:          IF (NITEMS.GT.1) CALL READF(SYMTOL1)
6714:          IF (NITEMS.GT.2) CALL READF(SYMTOL2)6710:          IF (NITEMS.GT.2) CALL READF(SYMTOL2)
6715:          IF (NITEMS.GT.3) CALL READF(SYMTOL3)6711:          IF (NITEMS.GT.3) CALL READF(SYMTOL3)
6716: 6712: 
6717:       ELSE IF (WORD .EQ. 'PRINT_MINDATA') THEN 
6718:          PRINT_MINDATA = .TRUE. 
6719:  
6720: !  Keyword and parameters for symmetrisation.6713: !  Keyword and parameters for symmetrisation.
6721: !6714: !
6722:       ELSE IF (WORD.EQ.'SYMMETRISE') THEN6715:       ELSE IF (WORD.EQ.'SYMMETRISE') THEN
6723:          SYMMETRIZE=.TRUE.6716:          SYMMETRIZE=.TRUE.
6724:          NCORE=06717:          NCORE=0
6725:          IF (NITEMS.GT.1) CALL READI(NSYMINTERVAL)6718:          IF (NITEMS.GT.1) CALL READI(NSYMINTERVAL)
6726:          IF (NITEMS.GT.2) CALL READF(SYMTOL1)6719:          IF (NITEMS.GT.2) CALL READF(SYMTOL1)
6727:          IF (NITEMS.GT.3) CALL READF(SYMTOL2)6720:          IF (NITEMS.GT.3) CALL READF(SYMTOL2)
6728:          IF (NITEMS.GT.4) CALL READF(SYMTOL3)6721:          IF (NITEMS.GT.4) CALL READF(SYMTOL3)
6729:          IF (NITEMS.GT.5) CALL READF(SYMTOL4)6722:          IF (NITEMS.GT.5) CALL READF(SYMTOL4)


r29933/mie_field.f90 2016-03-16 18:33:40.475131930 +0000 r29932/mie_field.f90 2016-03-16 18:33:41.911146696 +0000
  1: SUBROUTINE MIEF_INI()  1: SUBROUTINE MIEF_INI()
  2:   !  2:   !
  3:   USE COMMONS, ONLY : MIEF_FILENAME, MIEF_EPS, MIEF_SIG, &  3:   USE COMMONS, ONLY : MIEF_FILENAME, MIEF_EPS, MIEF_SIG, &
  4:        MIEF_NSITES, MIEF_N, MIEF_M, NSPECIES, MIEF_SITES, &  4:        MIEF_NSITES, MIEF_N, MIEF_M, NSPECIES, MIEF_SITES, &
  5:        MIEF_CUTT, MIEF_RCUT, MIEF_U_RCUT, MIEF_DUDR_RCUT, MYUNIT  5:        MIEF_CUTT, MIEF_RCUT, MIEF_U_RCUT, MIEF_DUDR_RCUT, MYUNIT
  6:   !  6:   !
  7:   IMPLICIT NONE  7:   IMPLICIT NONE
  8:   !  8:   !
  9:   LOGICAL :: YESNO  9:   LOGICAL :: YESNO
 10:   INTEGER :: I, J, K, L, EOF, LUNIT, GETUNIT 10:   INTEGER :: I, J, EOF, LUNIT, GETUNIT
 11:   DOUBLE PRECISION :: PREF, DUMMY, REP_TERM, ATT_TERM 11:   DOUBLE PRECISION :: PREF, DUMMY, REP_TERM, ATT_TERM
 12:   ! 12:   !
 13:   ALLOCATE(MIEF_SIG(NSPECIES(0))) 13:   ALLOCATE(MIEF_SIG(NSPECIES(0)))
 14:   ALLOCATE(MIEF_EPS(NSPECIES(0))) 14:   ALLOCATE(MIEF_EPS(NSPECIES(0)))
 15:   ! 15:   !
 16:   ! Process the file with Mie parameters and site coordinates 16:   ! Process the file with Mie parameters and site coordinates
 17:   ! 17:   !
 18:   INQUIRE(FILE=MIEF_FILENAME, EXIST=YESNO) 18:   INQUIRE(FILE=MIEF_FILENAME, EXIST=YESNO)
 19:   ! 19:   !
 20:   IF(YESNO) THEN 20:   IF(YESNO) THEN
 74:         MIEF_U_RCUT(I)=PREF*MIEF_EPS(I)*(REP_TERM-ATT_TERM) 74:         MIEF_U_RCUT(I)=PREF*MIEF_EPS(I)*(REP_TERM-ATT_TERM)
 75:         MIEF_DUDR_RCUT(I)=PREF*MIEF_EPS(I)*( -DBLE(MIEF_N)*REP_TERM + & 75:         MIEF_DUDR_RCUT(I)=PREF*MIEF_EPS(I)*( -DBLE(MIEF_N)*REP_TERM + &
 76:              DBLE(MIEF_M)*ATT_TERM ) / MIEF_RCUT 76:              DBLE(MIEF_M)*ATT_TERM ) / MIEF_RCUT
 77:      ENDDO 77:      ENDDO
 78:   ENDIF 78:   ENDIF
 79:   ! 79:   !
 80:   RETURN 80:   RETURN
 81:   ! 81:   !
 82: END SUBROUTINE MIEF_INI 82: END SUBROUTINE MIEF_INI
 83: ! 83: !
 84: SUBROUTINE MIEF(X,GRAD,EREAL,GRADT,HESST,STRESST) 84: SUBROUTINE MIEF(X,GRAD,EREAL,GRADT, STRESST)
 85:   ! 85:   !
 86:   USE COMMONS, ONLY : MIEF_N, MIEF_M, MIEF_NSITES, MIEF_SITES, & 86:   USE COMMONS, ONLY : MIEF_N, MIEF_M, MIEF_NSITES, MIEF_SITES, &
 87:        MIEF_EPS, MIEF_SIG, ATOMLISTS, INVATOMLISTS, VT, NATOMS, & 87:        MIEF_EPS, MIEF_SIG, ATOMLISTS, INVATOMLISTS, VT, NATOMS, &
 88:        NSPECIES, MIEF_PBCT, MIEF_BOX, MIEF_CUTT, MIEF_RCUT, & 88:        NSPECIES, MIEF_PBCT, MIEF_BOX, MIEF_CUTT, MIEF_RCUT, &
 89:        MIEF_U_RCUT, MIEF_DUDR_RCUT, STRESS 89:        MIEF_U_RCUT, MIEF_DUDR_RCUT, STRESS
 90:   USE MODHESS, ONLY : HESS 
 91:   ! 90:   !
 92:   IMPLICIT NONE 91:   IMPLICIT NONE
 93:   ! 92:   !
 94:   DOUBLE PRECISION, INTENT(IN) :: X(3*NATOMS) 93:   DOUBLE PRECISION, INTENT(IN) :: X(3*NATOMS)
 95:   DOUBLE PRECISION, INTENT(INOUT) :: GRAD(3*NATOMS) 94:   DOUBLE PRECISION, INTENT(INOUT) :: GRAD(3*NATOMS)
 96:   DOUBLE PRECISION, INTENT(INOUT) :: EREAL 95:   DOUBLE PRECISION, INTENT(INOUT) :: EREAL
 97:   LOGICAL, INTENT(IN) :: GRADT, HESST, STRESST 96:   LOGICAL, INTENT(IN) :: GRADT, STRESST
 98:   ! 97:   !
 99:   INTEGER :: I, J, K, L, I3, ITYPE 98:   INTEGER :: I, J, K, L, I3, ITYPE
100:   DOUBLE PRECISION :: PREF, DX(3),DIST2,IDIST,SIG,EPS,& 99:   DOUBLE PRECISION :: PREF, DX(3),DIST2,IDIST,SIG,EPS,&
101:        DUMMY, ATT_TERM, REP_TERM, DIST, DUDR, D2UDR2100:        DUMMY, ATT_TERM, REP_TERM, DIST
102:   !101:   !
103:   PREF = DBLE(MIEF_N)/DBLE(MIEF_N-MIEF_M)*&102:   PREF = DBLE(MIEF_N)/DBLE(MIEF_N-MIEF_M)*&
104:        (DBLE(MIEF_N)/DBLE(MIEF_M))**&103:        (DBLE(MIEF_N)/DBLE(MIEF_M))**&
105:        (DBLE(MIEF_M)/DBLE(MIEF_N-MIEF_M))104:        (DBLE(MIEF_M)/DBLE(MIEF_N-MIEF_M))
106:   !105:   !
107:   DO I=1,NATOMS ! Loop over ALL atoms106:   DO I=1,NATOMS ! Loop over ALL atoms
108:      !107:      !
109:      I3=3*(I-1)108:      I3=3*(I-1)
110:      ITYPE=INVATOMLISTS(I,1)109:      ITYPE=INVATOMLISTS(I,1)
111:      EPS=MIEF_EPS(ITYPE)*PREF110:      EPS=MIEF_EPS(ITYPE)*PREF
131:            DUMMY = EPS*(REP_TERM-ATT_TERM)130:            DUMMY = EPS*(REP_TERM-ATT_TERM)
132:            IF(MIEF_CUTT) THEN131:            IF(MIEF_CUTT) THEN
133:               DUMMY = DUMMY - MIEF_U_RCUT(ITYPE) - &132:               DUMMY = DUMMY - MIEF_U_RCUT(ITYPE) - &
134:                    (DIST-MIEF_RCUT)*MIEF_DUDR_RCUT(ITYPE)133:                    (DIST-MIEF_RCUT)*MIEF_DUDR_RCUT(ITYPE)
135:            ENDIF134:            ENDIF
136:            VT(I) = VT(I) + DUMMY135:            VT(I) = VT(I) + DUMMY
137:            EREAL = EREAL + DUMMY136:            EREAL = EREAL + DUMMY
138:            !137:            !
139:            IF(GRADT) THEN138:            IF(GRADT) THEN
140:               ! Compute 1/R*dU/dR139:               ! Compute 1/R*dU/dR
141:               DUDR = EPS*( -DBLE(MIEF_N)*REP_TERM + &140:               DUMMY = EPS*( -DBLE(MIEF_N)*REP_TERM + &
142:                    DBLE(MIEF_M)*ATT_TERM ) / DIST2141:                    DBLE(MIEF_M)*ATT_TERM ) / DIST2
143:               IF(MIEF_CUTT) DUDR = DUDR - MIEF_DUDR_RCUT(ITYPE)/DIST142:               IF(MIEF_CUTT) DUMMY = DUMMY - MIEF_DUDR_RCUT(ITYPE)/DIST
144:               DO K=1,3143:               DO K=1,3
145:                  GRAD(I3+K) = GRAD(I3+K) + DUDR*DX(K)144:                  GRAD(I3+K) = GRAD(I3+K) + DUMMY*DX(K)
146:                  IF(STRESST) THEN ! Accumulate local stresses145:                  IF(STRESST) THEN ! Accumulate local stresses
147:                     DO L=K,3146:                     DO L=K,3
148:                        STRESS(I,K,L) = STRESS(I,K,L) + &147:                        STRESS(I,K,L) = STRESS(I,K,L) + &
149:                             DUDR*DX(K)*DX(L)148:                             DUMMY*DX(K)*DX(L)
150:                     ENDDO149:                     ENDDO
151:                  ENDIF150:                  ENDIF
152:               ENDDO151:               ENDDO
153:               ! 
154:               IF(HESST) THEN 
155:                  ! Compute d2U/dR2                                      
156:                  D2UDR2 = EPS*( DBLE(MIEF_N*(MIEF_N+1))*REP_TERM - & 
157:                       DBLE(MIEF_M*(MIEF_M+1))*ATT_TERM ) / DIST2 
158:                  ! What we need is (D2UDR - DUDR/DIST) / DIST2 
159:                  D2UDR2 = (D2UDR2 - DUDR) / DIST2 
160:                  ! 
161:                  DO K=1,3 ! 1st coordinate of atom I 
162:                     DO L=K,3 ! 2nd coordinate of atom I 
163:                        DUMMY = DX(K)*DX(L) ! Reset DUMMY 
164:                        DUMMY = DUMMY*D2UDR2 
165:                        IF(K==L) DUMMY = DUMMY + DUDR 
166:                        ! Accumulate diagonal 3-blocks only 
167:                        HESS(I3+K, I3+L) = HESS(I3+K, I3+L) + DUMMY 
168:                     ENDDO 
169:                  ENDDO 
170:                  ! 
171:               ENDIF 
172:               ! 
173:            ENDIF152:            ENDIF
174:            !153:            !
175:         ENDIF ! Cutoff check154:         ENDIF ! Cutoff check
176:         !155:         !
177:      ENDDO ! Spanned sites156:      ENDDO
178:      !157:   ENDDO
179:      IF(HESST) THEN 
180:         ! Use symmetry to fill skipped entries in diagonal blocks 
181:         DO K=1,3 
182:            DO L=K+1,3 
183:               HESS(I3+L,I3+K) = HESS(I3+K, I3+L) 
184:            ENDDO 
185:         ENDDO 
186:      END IF 
187:      ! 
188:   ENDDO ! Spanned atoms 
189:   !158:   !
190:   IF(GRADT) THEN159:   IF(GRADT) THEN
191:      DO I=1,NSPECIES(0) ! Atom types160:      DO I=1,NSPECIES(0) ! Atom types
192:         DO J=1,ATOMLISTS(I,2,0) ! group 2 only161:         DO J=1,ATOMLISTS(I,2,0) ! group 2 only
193:            K=3*(ATOMLISTS(I,2,J)-1)162:            K=3*(ATOMLISTS(I,2,J)-1)
194:            GRAD(K+1) = 0.0D0 ! Reset to zero163:            GRAD(K+1) = 0.0D0 ! Reset to zero
195:            GRAD(K+2) = 0.0D0164:            GRAD(K+2) = 0.0D0
196:            GRAD(K+3) = 0.0D0165:            GRAD(K+3) = 0.0D0
197:         ENDDO166:         ENDDO
198:      END DO167:      END DO


r29933/MLJ.f90 2016-03-16 18:33:39.687123826 +0000 r29932/MLJ.f90 2016-03-16 18:33:41.115138512 +0000
 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,ATOMLISTS, NSPECIES, STRESS 29:   USE COMMONS, ONLY : NATOMS,VT,ATOMLISTS, NSPECIES, STRESS
 30:   USE MODHESS, ONLY : HESS 
 31:   USE POT_PARAMS, ONLY : MLJ_EPS, MLJ_SIG 30:   USE POT_PARAMS, ONLY : MLJ_EPS, MLJ_SIG
 32:   ! 31:   !
 33:   IMPLICIT NONE 32:   IMPLICIT NONE
 34:   ! 33:   !
 35:   LOGICAL, INTENT (IN) :: GRADT,HESST,STRESST 34:   LOGICAL, INTENT (IN) :: GRADT,HESST,STRESST
 36:   DOUBLE PRECISION, INTENT (IN) :: X(3*NATOMS) 35:   DOUBLE PRECISION, INTENT (IN) :: X(3*NATOMS)
 37:   DOUBLE PRECISION, INTENT (INOUT) :: POT, GRAD(3*NATOMS) 36:   DOUBLE PRECISION, INTENT (INOUT) :: POT, GRAD(3*NATOMS)
 38:   ! 37:   !
 39:   LOGICAL :: SELF 38:   LOGICAL :: SELF
 40:   INTEGER :: I,J,J1,J13,J2,J23,T1,T2,G1,G2,G2START,& 39:   INTEGER :: I,J,J1,J13,J2,J23,T1,T2,G1,G2,G2START,&
 41:        K1,K2,LI1,LI2,UI1,UI2 40:        K1,K2,LI1,LI2,UI1,UI2
 42:   DOUBLE PRECISION :: DIST2, DX(3), IR6, IR12, DUMMY, & 41:   DOUBLE PRECISION :: DIST2, DX(3), IR6, IR12, DUMMY, &
 43:        SIG2, EPS4, EPS24, IGRAD(3), DUDR, D2UDR2 42:        SIG2, EPS4, EPS24, IGRAD(3)
 44:   ! 43:   !
 45:   ! Zero the potential and the gradient 44:   ! Zero the potential and the gradient
 46:   VT(1:NATOMS) = 0.0D0 45:   VT(1:NATOMS) = 0.0D0
 47:   POT = 0.0D0 46:   POT = 0.0D0
 48:   IF(GRADT) GRAD(:) = 0.0D0 47:   IF(GRADT) GRAD(:) = 0.0D0
 49:   IF(HESST) HESS(:,:) = 0.0D0 
 50:   IF(STRESST) STRESS(:,:,:) = 0.0D0 48:   IF(STRESST) STRESS(:,:,:) = 0.0D0
 51:   ! 49:   !
 52:   ! Double loop over atom types 50:   ! Double loop over atom types
 53:   LT1: DO T1=1,NSPECIES(0) 51:   LT1: DO T1=1,NSPECIES(0)
 54:      LT2: DO T2=1,T1  52:      LT2: DO T2=1,T1 
 55:         ! 53:         !
 56:         ! This can be precomputed elsewhere... 54:         ! This can be precomputed elsewhere...
 57:         SIG2 = MLJ_SIG(T1,T2)*MLJ_SIG(T1,T2) 55:         SIG2 = MLJ_SIG(T1,T2)*MLJ_SIG(T1,T2)
 58:         EPS4 = 4.0D0*MLJ_EPS(T1,T2) 56:         EPS4 = 4.0D0*MLJ_EPS(T1,T2)
 59:         IF(GRADT) EPS24 = 24.0D0*MLJ_EPS(T1,T2) 57:         IF(GRADT) EPS24 = 24.0D0*MLJ_EPS(T1,T2)
110:                     IR12 = IR6*IR6108:                     IR12 = IR6*IR6
111:                     DUMMY = EPS4*(IR12 - IR6)109:                     DUMMY = EPS4*(IR12 - IR6)
112:                     !WRITE(*,*) "DUMMY FOR POT", DUMMY110:                     !WRITE(*,*) "DUMMY FOR POT", DUMMY
113:                     ! -----------------------------------------111:                     ! -----------------------------------------
114:                     !112:                     !
115:                     VT(J1) = VT(J1) + DUMMY113:                     VT(J1) = VT(J1) + DUMMY
116:                     VT(J2) = VT(J2) + DUMMY114:                     VT(J2) = VT(J2) + DUMMY
117:                     POT = POT + DUMMY115:                     POT = POT + DUMMY
118:                     !116:                     !
119:                     IF(GRADT) THEN ! Calculate gradient117:                     IF(GRADT) THEN ! Calculate gradient
120:                        ! Compute 1/R*dUdR118:                        !
121:                        DUDR = EPS24*(IR6 - 2.0D0*IR12)/DIST2119:                        DUMMY = EPS24*(IR6 - 2.0D0*IR12)/DIST2
122:                        !WRITE(*,*) "DUMMY FOR GRAD", DUMMY120:                        !WRITE(*,*) "DUMMY FOR GRAD", DUMMY
123:                        DO I = 1,3121:                        DO I = 1,3
124:                           IGRAD(I) = DUDR*DX(I)122:                           IGRAD(I) = DUMMY*DX(I)
125:                           GRAD(J13+I) = GRAD(J13+I) + IGRAD(I)123:                           GRAD(J13+I) = GRAD(J13+I) + IGRAD(I)
126:                           GRAD(J23+I) = GRAD(J23+I) - IGRAD(I)124:                           GRAD(J23+I) = GRAD(J23+I) - IGRAD(I)
127:                        ENDDO125:                        ENDDO
128:                        !126:                        !
129:                        IF(STRESST) THEN127:                        IF(STRESST) THEN
130:                           DO I=1,3128:                           DO I=1,3
131:                              DO J=I,3129:                              DO J=I,3
132:                                 DUMMY = IGRAD(I)*DX(J)130:                                 DUMMY = IGRAD(I)*DX(J)
133:                                 STRESS(J1,I,J) = STRESS(J1,I,J) + &131:                                 STRESS(J1,I,J) = STRESS(J1,I,J) + &
134:                                      DUMMY132:                                      DUMMY
135:                                 STRESS(J2,I,J) = STRESS(J2,I,J) + &133:                                 STRESS(J2,I,J) = STRESS(J2,I,J) + &
136:                                      DUMMY134:                                      DUMMY
137:                                 STRESS(0,I,J) = STRESS(0,I,J) + &135:                                 STRESS(0,I,J) = STRESS(0,I,J) + &
138:                                      DUMMY136:                                      DUMMY
139:                              ENDDO137:                              ENDDO
140:                           ENDDO138:                           ENDDO
141:                        ENDIF139:                        ENDIF
142:                        !140:                        !
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+I, J13+J) = 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:                        ! 
186:                     ENDIF141:                     ENDIF
187:                     !142:                     !
188:                  ENDDO LA2 ! Double loop over atoms143:                  ENDDO LA2 ! Double loop over atoms
189:               ENDDO LA1144:               ENDDO LA1
190:               !145:               !
191:            ENDDO LG2 ! Double loop over groups146:            ENDDO LG2 ! Double loop over groups
192:         ENDDO LG1147:         ENDDO LG1
193:         !148:         !
194:      ENDDO LT2 ! Double loop over types149:      ENDDO LT2 ! Double loop over types
195:   ENDDO LT1150:   ENDDO LT1
212:      !            T1, G1, ATOMLISTS(T1,G1,0)167:      !            T1, G1, ATOMLISTS(T1,G1,0)
213:      !       DO J2=1,ATOMLISTS(T1,G1,0)168:      !       DO J2=1,ATOMLISTS(T1,G1,0)
214:      !          J1=ATOMLISTS(T1,G1,J2)169:      !          J1=ATOMLISTS(T1,G1,J2)
215:      !          J13=3*(J1-1)170:      !          J13=3*(J1-1)
216:      !          WRITE(*,'(A,1X, I4, 4(1X,F12.6))') "BLJ_CLUST>", &171:      !          WRITE(*,'(A,1X, I4, 4(1X,F12.6))') "BLJ_CLUST>", &
217:      !               J1, VT(J1), GRAD(J13+1), GRAD(J13+2), GRAD(J13+3)172:      !               J1, VT(J1), GRAD(J13+1), GRAD(J13+2), GRAD(J13+3)
218:      !       ENDDO173:      !       ENDDO
219:      !    ENDDO174:      !    ENDDO
220:      ! ENDDO175:      ! ENDDO
221:      !176:      !
222:      IF(STRESST) THEN177:   ENDIF
223:         DO I=1,3178:   !
224:            DO J=I,3179:   IF(STRESST) THEN
225:               STRESS(0,J,I) = STRESS(0,I,J) ! Impose symmetry180:      DO I=1,3
226:            ENDDO181:         DO J=I,3
227:         ENDDO182:            STRESS(0,J,I) = STRESS(0,I,J) ! Impose symmetry
228:      ENDIF 
229:      ! 
230:      IF(HESST) THEN 
231:         ! Use symmetry to fill the skipped entries in diagonal blocks 
232:         DO J1=1,NATOMS 
233:            DO I = 1,3 ! coordinates of atom 1 
234:               DO J = I+1,3 ! coordinates of atom 2  
235:                  HESS(J1+J, J1+I) = HESS(J1+I, J1+J) 
236:               ENDDO 
237:            ENDDO 
238:         ENDDO183:         ENDDO
239:         !184:      ENDDO
240:      ENDIF 
241:      ! 
242:   ENDIF185:   ENDIF
243:   !186:   !
244:   RETURN187:   RETURN
245:   !188:   !
246: END SUBROUTINE MLJ189: END SUBROUTINE MLJ


r29933/potential.f90 2016-03-16 18:33:40.687134110 +0000 r29932/potential.f90 2016-03-16 18:33:42.107148712 +0000
1234: 1234: 
1235:    IF (K_COMP == 0.0D0) COMPON=.FALSE.1235:    IF (K_COMP == 0.0D0) COMPON=.FALSE.
1236:    IF (COMPON) CALL COMPRESS(X, GRAD, EREAL, GRADT)1236:    IF (COMPON) CALL COMPRESS(X, GRAD, EREAL, GRADT)
1237: 1237: 
1238: ! Apply twist forces, if appropriate1238: ! Apply twist forces, if appropriate
1239:    IF (TWISTT) THEN1239:    IF (TWISTT) THEN
1240:       CALL TWIST(X(1:3*NATOMS), NATOMS, GRAD(1:3*NATOMS), EREAL, GRADT)1240:       CALL TWIST(X(1:3*NATOMS), NATOMS, GRAD(1:3*NATOMS), EREAL, GRADT)
1241:    END IF1241:    END IF
1242: 1242: 
1243: ! ds656> Apply a substrate field for keyword MIE_FIELD1243: ! ds656> Apply a substrate field for keyword MIE_FIELD
1244:    IF(MIEFT) CALL MIEF(X,GRAD,EREAL,GRADT,SECT,.FALSE.)1244:    IF(MIEFT) CALL MIEF(X,GRAD,EREAL,GRADT,.FALSE.)
1245: 1245: 
1246: !1246: !
1247: ! js850> apply harmonic field for keyword HARMONICF1247: ! js850> apply harmonic field for keyword HARMONICF
1248: !1248: !
1249:    IF ( HARMONICF ) THEN1249:    IF ( HARMONICF ) THEN
1250:       DO J1=1, NATOMS1250:       DO J1=1, NATOMS
1251:          IF (HARMONICFLIST(J1)) THEN1251:          IF (HARMONICFLIST(J1)) THEN
1252: !            WRITE(*,'(1x, 6F25.16)') X(1), X(2), X(3), X0(1), X0(2), X0(3)1252: !            WRITE(*,'(1x, 6F25.16)') X(1), X(2), X(3), X0(1), X0(2), X0(3)
1253: !            WRITE(*,*) "#V(1) ", GRAD((J1-1)*3+1)1253: !            WRITE(*,*) "#V(1) ", GRAD((J1-1)*3+1)
1254:             CALL HARMONICFIELD( X(3*J1-2), HARMONICR0(3*J1-2), GRAD(3*J1-2), EREAL)1254:             CALL HARMONICFIELD( X(3*J1-2), HARMONICR0(3*J1-2), GRAD(3*J1-2), EREAL)


r29933/stress.f90 2016-03-16 18:33:40.887136166 +0000 r29932/stress.f90 2016-03-16 18:33:42.299150685 +0000
  6:   IMPLICIT NONE  6:   IMPLICIT NONE
  7:   !  7:   !
  8:   DOUBLE PRECISION, INTENT(IN) :: X(3*NATOMS)  8:   DOUBLE PRECISION, INTENT(IN) :: X(3*NATOMS)
  9:   DOUBLE PRECISION, INTENT(OUT) :: EPOT  9:   DOUBLE PRECISION, INTENT(OUT) :: EPOT
 10:   ! 10:   !
 11:   INTEGER :: I,J 11:   INTEGER :: I,J
 12:   DOUBLE PRECISION :: GRAD(3*NATOMS) 12:   DOUBLE PRECISION :: GRAD(3*NATOMS)
 13:   ! 13:   !
 14:   IF(.NOT.STRESST) THEN 14:   IF(.NOT.STRESST) THEN
 15:      WRITE(MYUNIT,'(A)') & 15:      WRITE(MYUNIT,'(A)') &
 16:           'calc_stress> WARNING: Should not be here!' 16:           'calc_stress> Should not be here!'
 17:      !RETURN 17:      RETURN
 18:   ELSE 18:   ELSE
 19:      STRESS(:,:,:) = 0.0D0 ! Initialise 19:      STRESS(:,:,:) = 0.0D0 ! Initialise
 20:   ENDIF 20:   ENDIF
 21:   ! 21:   !
 22:   IF(MGUPTAT) THEN 22:   IF(MGUPTAT) THEN
 23:      CALL MGUPTA(X, GRAD, EPOT, .TRUE., .FALSE., .TRUE.) 23:      CALL MGUPTA(X, GRAD, EPOT, .TRUE., .FALSE., .TRUE.)
 24:   ELSEIF(MLJT) THEN 24:   ELSEIF(MLJT) THEN
 25:      CALL MLJ(X, GRAD, EPOT, .TRUE., .FALSE., .TRUE.) 25:      CALL MLJ(X, GRAD, EPOT, .TRUE., .FALSE., .TRUE.)
 26:   ELSE 26:   ELSE
 27:      WRITE(MYUNIT,'(A)') & 27:      WRITE(MYUNIT,'(A)') &
 28:           'calc_stress> Stress calculation not implemented for current potential.' 28:           'calc_stress> Stress calculation not implemented for current potential.'
 29:   ENDIF 29:   ENDIF
 30:   ! 30:   !
 31:   IF(MIEFT) CALL MIEF(X,GRAD,EPOT,.TRUE.,.FALSE.,.TRUE.) 31:   IF(MIEFT) CALL MIEF(X,GRAD,EPOT,.TRUE.,.TRUE.)
 32:   ! 32:   !
 33:   WRITE(MYUNIT,'(A)') 'stress> Overall stress tensor:' 33:   WRITE(MYUNIT,'(A)') 'stress> Overall stress tensor:'
 34:   DO I=1,3 34:   DO I=1,3
 35:      WRITE(MYUNIT,'(3(1X,E15.8))') (STRESS(0,I,J),J=1,3) 35:      WRITE(MYUNIT,'(3(1X,E15.8))') (STRESS(0,I,J),J=1,3)
 36:   ENDDO 36:   ENDDO
 37:   ! 37:   !
 38:   RETURN 38:   RETURN
 39:   ! 39:   !
 40: END SUBROUTINE CALC_STRESS 40: END SUBROUTINE CALC_STRESS


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0