hdiff output

r31913/commons.f90 2017-02-15 08:30:16.480718991 +0000 r31912/commons.f90 2017-02-15 08:30:17.648734603 +0000
112:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, POLYT, SANDBOXT, &112:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, POLYT, SANDBOXT, &
113:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &113:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
114:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &114:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
115:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &115:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
116:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &116:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &
117:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 117:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 
118:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &118:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &
119:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &119:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
120:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &120:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &
121:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &121:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &
122:      &        DUMPMQT, MLQT, MLQPROB, LJADD2T, MLPVB3T, NOREGBIAS, PYADDT, PYADD2T122:      &        DUMPMQT, MLQT, MLQPROB, LJADD2T, MLPVB3T, NOREGBIAS, PYADDT
123: !123: !
124:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 124:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 
125:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)125:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)
126:       DOUBLE PRECISION, ALLOCATABLE :: SPECMASS(:) 126:       DOUBLE PRECISION, ALLOCATABLE :: SPECMASS(:) 
127: 127: 
128: ! csw34> FREEZEGROUP variables128: ! csw34> FREEZEGROUP variables
129: !129: !
130:       INTEGER :: GROUPCENTRE130:       INTEGER :: GROUPCENTRE
131:       DOUBLE PRECISION :: GROUPRADIUS131:       DOUBLE PRECISION :: GROUPRADIUS
132:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE132:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE
620:       DOUBLE PRECISION, ALLOCATABLE :: EPSALPHA(:), DISTREF(:)620:       DOUBLE PRECISION, ALLOCATABLE :: EPSALPHA(:), DISTREF(:)
621:       DOUBLE PRECISION, ALLOCATABLE :: CPDISTREF(:)621:       DOUBLE PRECISION, ALLOCATABLE :: CPDISTREF(:)
622:       DOUBLE PRECISION, ALLOCATABLE :: MIN1REDO(:), MIN2REDO(:)622:       DOUBLE PRECISION, ALLOCATABLE :: MIN1REDO(:), MIN2REDO(:)
623:       DOUBLE PRECISION, ALLOCATABLE :: CONDISTREF(:), CONDISTREFLOCAL(:), CONDISTREFON(:), CONDISTREFLOCALON(:)623:       DOUBLE PRECISION, ALLOCATABLE :: CONDISTREF(:), CONDISTREFLOCAL(:), CONDISTREFON(:), CONDISTREFLOCALON(:)
624:       DOUBLE PRECISION, ALLOCATABLE :: CONCUT(:), CONCUTLOCAL(:)624:       DOUBLE PRECISION, ALLOCATABLE :: CONCUT(:), CONCUTLOCAL(:)
625:       DOUBLE PRECISION, ALLOCATABLE :: CONDISTREFFIX(:), CONCUTFIX(:)625:       DOUBLE PRECISION, ALLOCATABLE :: CONDISTREFFIX(:), CONCUTFIX(:)
626:       DOUBLE PRECISION, ALLOCATABLE :: CPCONDISTREF(:)626:       DOUBLE PRECISION, ALLOCATABLE :: CPCONDISTREF(:)
627:       DOUBLE PRECISION, ALLOCATABLE :: CONGEOM(:,:)627:       DOUBLE PRECISION, ALLOCATABLE :: CONGEOM(:,:)
628:       DOUBLE PRECISION, ALLOCATABLE :: LJADDEPS(:,:)628:       DOUBLE PRECISION, ALLOCATABLE :: LJADDEPS(:,:)
629:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)629:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)
630:       DOUBLE PRECISION, ALLOCATABLE :: PYADDREP(:,:), PYADDATT(:,:) 
631:       LOGICAL, ALLOCATABLE :: CONACTIVE(:)630:       LOGICAL, ALLOCATABLE :: CONACTIVE(:)
632:       LOGICAL, ALLOCATABLE :: ATOMACTIVE(:)631:       LOGICAL, ALLOCATABLE :: ATOMACTIVE(:)
633:       INTEGER, ALLOCATABLE :: ORDERI(:), ORDERJ(:), REPPOW(:)632:       INTEGER, ALLOCATABLE :: ORDERI(:), ORDERJ(:), REPPOW(:)
634:       INTEGER, ALLOCATABLE :: CONI(:), CONJ(:), CONION(:), CONJON(:)633:       INTEGER, ALLOCATABLE :: CONI(:), CONJ(:), CONION(:), CONJON(:)
635:       INTEGER, ALLOCATABLE :: CONIFIX(:), CONJFIX(:), REPIFIX(:), REPJFIX(:)634:       INTEGER, ALLOCATABLE :: CONIFIX(:), CONJFIX(:), REPIFIX(:), REPJFIX(:)
636:       INTEGER, ALLOCATABLE :: REPI(:), REPJ(:)635:       INTEGER, ALLOCATABLE :: REPI(:), REPJ(:)
637:       INTEGER, ALLOCATABLE :: CPCONI(:), CPCONJ(:)636:       INTEGER, ALLOCATABLE :: CPCONI(:), CPCONJ(:)
638:       INTEGER, ALLOCATABLE :: CPREPI(:), CPREPJ(:)637:       INTEGER, ALLOCATABLE :: CPREPI(:), CPREPJ(:)
639:       DOUBLE PRECISION, ALLOCATABLE :: REPCUT(:), NREPCUT(:), CPREPCUT(:), REPCUTFIX(:)638:       DOUBLE PRECISION, ALLOCATABLE :: REPCUT(:), NREPCUT(:), CPREPCUT(:), REPCUTFIX(:)
640:       INTEGER, ALLOCATABLE :: NREPI(:), NREPJ(:)639:       INTEGER, ALLOCATABLE :: NREPI(:), NREPJ(:)


r31913/gay-berne.f90 2017-02-15 08:30:16.704721986 +0000 r31912/gay-berne.f90 2017-02-15 08:30:17.872737597 +0000
1046:                  DFP(2,3)=DF2PI31046:                  DFP(2,3)=DF2PI3
1047:                  DFP(2,4)=DF2PJ11047:                  DFP(2,4)=DF2PJ1
1048:                  DFP(2,5)=DF2PJ21048:                  DFP(2,5)=DF2PJ2
1049:                  DFP(2,6)=DF2PJ31049:                  DFP(2,6)=DF2PJ3
1050:                  DFDR(1,:)=DF1DR(:)1050:                  DFDR(1,:)=DF1DR(:)
1051:                  DFDR(2,:)=DF2DR(:)1051:                  DFDR(2,:)=DF2DR(:)
1052:                  LJ1(1)=RHO11052:                  LJ1(1)=RHO1
1053:                  LJ1(2)=RHO21053:                  LJ1(2)=RHO2
1054:  !     CALCULATE PY POTENTIAL ENERGY1054:  !     CALCULATE PY POTENTIAL ENERGY
1055: 123 CONTINUE 1055: 123 CONTINUE 
 1056:         VDUMMY = 0.0D0
 1057:                  
1056:                 VDUMMY=0.0D01058:                 VDUMMY=0.0D0
1057:                IF(LJSITE) THEN  1059:                IF(LJSITE) THEN  
1058:                 do k=1,maxinteractions ! k=1 -- interaction between repulsive primary 'apex' sites1060:                 do k=1,maxinteractions ! k=1 -- interaction between repulsive primary 'apex' sites
1059:                          ! k=2 and k=3 -- interaction between secondary and primary 'apex' sites1061:                          ! k=2 and k=3 -- interaction between secondary and primary 'apex' sites
1060:                          ! k=4 -- interaction between secondary 'apex' sites (normal LJ interaction) 1062:                          ! k=4 -- interaction between secondary 'apex' sites (normal LJ interaction) 
1061:                         ! trying to modify code to allow for binary systems. 1063:                         ! trying to modify code to allow for binary systems. 
1062:                         ! apex site heights will be defined in absolute units,1064:                         ! apex site heights will be defined in absolute units,
1063:                         ! hence PYA1bin(J1,1) etc. will be removed from below1065:                         ! hence PYA1bin(J1,1) etc. will be removed from below
1064: 1066: 
1065:                    IF(k==1) THEN1067:                    IF(k==1) THEN
1132:                VT(J1) = VT(J1) + VDUMMY1134:                VT(J1) = VT(J1) + VDUMMY
1133:                VT(J2) = VT(J2) + VDUMMY        ! pair potentials1135:                VT(J2) = VT(J2) + VDUMMY        ! pair potentials
1134: 1136: 
1135:             dVDUMMY(:) = 0.0D01137:             dVDUMMY(:) = 0.0D0
1136: 1138: 
1137: !     CALCULATE GRADIENT1139: !     CALCULATE GRADIENT
1138: 1140: 
1139:              FIJ = 0.0D01141:              FIJ = 0.0D0
1140:              TIJ = 0.0D01142:              TIJ = 0.0D0
1141:              TJI = 0.0D01143:              TJI = 0.0D0
 1144: 
 1145:              FIJ = 0.0D0
 1146:              TIJ = 0.0D0
 1147:              TJI = 0.0D0
1142:              dVDUMMY(:)=0.0D01148:              dVDUMMY(:)=0.0D0
1143: 1149: 
1144:            IF(LJSITE) THEN1150:            IF(LJSITE) THEN
1145:             do k=1,maxinteractions1151:             do k=1,maxinteractions
1146:              do j=1,31152:              do j=1,3
1147:                dVDUMMY(j) = dVDUMMY(j) + 4.0D0*epsilon1(k,J1,J2)*(12.0D0*LLJ(11,k)*dLLJ1(k,j))*term2(k) + &1153:                dVDUMMY(j) = dVDUMMY(j) + 4.0D0*epsilon1(k,J1,J2)*(12.0D0*LLJ(11,k)*dLLJ1(k,j))*term2(k) + &
1148:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(12,k)*term3(k)*drlj(k,j)! extra LJ site derivatives, currently1154:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(12,k)*term3(k)*drlj(k,j)! extra LJ site derivatives, currently
1149:                                                                                 ! only repulsive (LJ12)1155:                                                                                 ! only repulsive (LJ12)
1150:                dVDUMMY(j) = dVDUMMY(j) - attr(k)*(4.0D0*epsilon1(k,J1,J2)*(6.0D0*LLJ(5,k)*dLLJ1(k,j))*term2(k) + &1156:                dVDUMMY(j) = dVDUMMY(j) - attr(k)*(4.0D0*epsilon1(k,J1,J2)*(6.0D0*LLJ(5,k)*dLLJ1(k,j))*term2(k) + &
1151:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(6,k)*term3(k)*drlj(k,j)) ! attractive secondary apex site1157:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(6,k)*term3(k)*drlj(k,j)) ! attractive secondary apex site
1200:                 IF(FROZEN(J1)) THEN1206:                 IF(FROZEN(J1)) THEN
1201:                         G(J2-2)=0.0D01207:                         G(J2-2)=0.0D0
1202:                         G(J2-1)=0.0D01208:                         G(J2-1)=0.0D0
1203:                         G(J2)=0.0D01209:                         G(J2)=0.0D0
1204:                 END IF1210:                 END IF
1205:         END DO1211:         END DO
1206:         RETURN1212:         RETURN
1207:       END SUBROUTINE PYGPERIODICADD1213:       END SUBROUTINE PYGPERIODICADD
1208: 1214: 
1209: 1215: 
1210:  
1211: ! PY potential, dc430's implementation 
1212: ! with PBC and continuous cutoff added 
1213: ! plus extra LJ site 
1214: ! 
1215: ! This version is addressable with interactions scaled by the input epsilon matrix PYADDEPS 
1216: ! 
1217:       SUBROUTINE PYGPERIODICADD2 (X, G, ENERGY, GTEST) 
1218:  
1219:        use commons, only: BOXLX,BOXLY,BOXLZ,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PCUTOFF,PARAMONOVCUTOFF,& 
1220:                 &         natoms,pysignot,pyepsnot,radift,LJSITE,&  
1221:                 &         MAXINTERACTIONS,VT,& 
1222:                 &         FROZEN, PYADDREP, PYADDATT, NADDTARGET 
1223: !                &         PSIGMAATTR, PEPSILONATTR, MYUNIT, PYBINARYTYPE1, PYBINARYT,& 
1224: !                &         PSCALEFAC2, PSCALEFAC1, PEPSILON1, BLJSITE, PYA2BIN, PYA1BIN 
1225:  
1226:        use pymodule 
1227:       IMPLICIT NONE 
1228:  
1229:       DOUBLE PRECISION :: X(3*NATOMS), G(3*NATOMS) 
1230:       DOUBLE PRECISION :: RI(3), RJ(3), RIJ(3), NR(3), RIJSQ, ABSRIJ, P(3) 
1231:       DOUBLE PRECISION :: AE1(3,3), BE1(3,3), AE2(3,3), BE2(3,3), APB(3,3), APBINV(3,3) 
1232:       DOUBLE PRECISION :: FCNT1, FCNT2, SRTFI1, SRTFI2, SRTFI(2), FMIN, LAMDAC1, LAMDAC2, ENERGY 
1233:       DOUBLE PRECISION :: RHO1, RHO1SQ, RHO16, RHO112, RHO2, RHO2SQ, RHO26 
1234:       DOUBLE PRECISION :: FCTR1, FCTR2, DVDF1, DVDF2  
1235:       DOUBLE PRECISION :: DF1PI1, DF1PI2, DF1PI3, DF2PI1, DF2PI2, DF2PI3 
1236:       DOUBLE PRECISION :: DF1PJ1, DF1PJ2, DF1PJ3, DF2PJ1, DF2PJ2, DF2PJ3  
1237:       DOUBLE PRECISION :: RMI(3,3), RMJ(3,3) 
1238:       DOUBLE PRECISION :: DPI1RM(3,3), DPI2RM(3,3), DPI3RM(3,3), DPJ1RM(3,3), DPJ2RM(3,3), DPJ3RM(3,3) 
1239:       DOUBLE PRECISION :: DF1DR(3), DF2DR(3), DG1DR(3), DG2DR(3) 
1240:       DOUBLE PRECISION :: ARIBRJ(3), XC(3), XCMRI(3), XCMRJ(3), FIJ(3), TIJ(3), TJI(3) 
1241:       double precision :: D1ABEZ(3,3), D2ABEZ(3,3), D3ABEZ(3,3), D1ABE(3,3), D2ABE(3,3), D3ABE(3,3)  
1242:       LOGICAL          :: GTEST 
1243: !     DOUBLE PRECISION :: E(3,3), ESQ(3,3), DE1(3,3), DE2(3,3), DE3(3,3), THETA, THETA2, CT, ST, RM(3,3) 
1244:  
1245: ! sf344 additions 
1246:  
1247:       DOUBLE PRECISION :: CLJ1(2),CLJ2(2),CLJ3(2),CLJ4(2),CLJ5(2),CLJ6(2),CLJ7(2),CLJ8(2),CLJ11(2),CLJ12(2), & 
1248:                           & CLJ13(2),CLJ14(2),DFDR(2,3),DFP(2,6),dr(6),dCLJ1(2,12),dVDUMMY(12),LJ1(2),DUMMY,DUMMY1,DUMMY2,& 
1249:                           & dLJ1(2,12),VDUMMY 
1250:       DOUBLE PRECISION :: term2(maxinteractions),term3(maxinteractions), & 
1251:                           & xlj(maxinteractions,2,3),rljvec(maxinteractions,3),rljunitvec(maxinteractions,3),& 
1252:                           & drlj(maxinteractions,12),rlj(maxinteractions),rlj2(maxinteractions) 
1253:       DOUBLE PRECISION :: LLJ(12,maxinteractions), dLLJ1(maxinteractions,12) 
1254:       INTEGER          :: k, MJ1, MJ2 
1255:  
1256:      VT(1:NATOMS/2)=0.0D0 
1257:   
1258:      term2(:)=1.0D0 
1259:      term3(:)=0.0D0 
1260:        
1261:          ENERGY = 0.D0 
1262:          G(:)   = 0.D0 
1263:          
1264:         DO J1=1, REALNATOMS 
1265:             J3      = 3*J1 
1266:             J5      = OFFSET + J3 
1267:             RI      = X(J3-2:J3) 
1268:             P       = X(J5-2:J5) 
1269:             angle=sqrt(dot_product(P,P)) 
1270:             if(angle>twopi) then 
1271: ! normalise angle-axis coordinates 
1272:                 X(J5-2:J5)=X(J5-2:J5)/angle 
1273:                 do 
1274:                   angle=angle-twopi 
1275:                   if(angle<2*pi) exit 
1276:                 end do 
1277: ! multiply with new angle 
1278:                 X(J5-2:J5)=X(J5-2:J5)*angle 
1279:             end if 
1280:  
1281:             CALL RMDRVT(P, RMIvec(J1,:,:), DPI1RMvec(J1,:,:), DPI2RMvec(J1,:,:), DPI3RMvec(J1,:,:), GTEST) 
1282:  
1283:         END DO 
1284:  
1285:          DO J1 = 1, REALNATOMS 
1286:            
1287:             MJ1=MOD(J1-1,NADDTARGET)+1 
1288:             J3      = 3*J1 
1289:             J5      = OFFSET + J3 
1290:             RI      = X(J3-2:J3) 
1291:             P       = X(J5-2:J5) 
1292: !     ROTATION MATRIX 
1293:  
1294: !            CALL RMDRVT(P, RMI, DPI1RM, DPI2RM, DPI3RM, GTEST) 
1295:             RMI(:,:)=RMIvec(J1,:,:) 
1296:             DPI1RM(:,:)=DPI1RMvec(J1,:,:) 
1297:             DPI2RM(:,:)=DPI2RMvec(J1,:,:) 
1298:             DPI3RM(:,:)=DPI3RMvec(J1,:,:) 
1299:  
1300:             AE1 = MATMUL(RMI,(MATMUL(AEZR1(J1,:,:),(TRANSPOSE(RMI))))) 
1301:  
1302:             IF (RADIFT) THEN 
1303:  
1304:                AE2 = MATMUL(RMI,(MATMUL(AEZR2(J1,:,:),(TRANSPOSE(RMI)))))          
1305:  
1306:             ENDIF 
1307:  
1308: !     BEGIN INNER LOOP OVER PARTICLES 
1309:  
1310:             DO J2 = J1 + 1, REALNATOMS 
1311:  
1312:                MJ2=MOD(J2-1,NADDTARGET)+1 
1313:  
1314:                J4     = 3*J2 
1315:                J6     = OFFSET + J4 
1316:                RJ     = X(J4-2:J4)  
1317:                P      = X(J6-2:J6) 
1318:  
1319: !     ROTATION MATRIX 
1320:  
1321: !               CALL RMDRVT(P, RMJ, DPJ1RM, DPJ2RM, DPJ3RM, GTEST)                
1322:                RMJ(:,:)=RMIvec(J2,:,:) 
1323:                DPJ1RM(:,:)=DPI1RMvec(J2,:,:) 
1324:                DPJ2RM(:,:)=DPI2RMvec(J2,:,:) 
1325:                DPJ3RM(:,:)=DPI3RMvec(J2,:,:) 
1326:       
1327:                BE1 = MATMUL(RMJ,(MATMUL(AEZR1(J2,:,:),(TRANSPOSE(RMJ))))) 
1328:  
1329:                IF (RADIFT) THEN 
1330:     
1331:                   BE2 = MATMUL(RMJ,(MATMUL(AEZR2(J2,:,:),(TRANSPOSE(RMJ))))) 
1332:  
1333:                ENDIF 
1334:  
1335: !     CALCULATE SEPARATION 
1336:  
1337:                RIJ    = RI - RJ 
1338:                RIJSQ  = DOT_PRODUCT(RIJ,RIJ) 
1339:                ABSRIJ = DSQRT(RIJSQ) 
1340:                NR     = RIJ / ABSRIJ 
1341:  
1342:                 IF(PARAMONOVCUTOFF.AND.RIJSQ>cut) GOTO 124 
1343:  
1344: !     CALCULATE ECF 
1345:  
1346:                CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE1, BE1, RIJ, LAMDAC1, FMIN) 
1347:  
1348:                FCNT1   = - FMIN 
1349:                SRTFI1  = 1.D0 / DSQRT(FCNT1) 
1350:                APB     = LAMDAC1 * AE1 + (1.D0 - LAMDAC1) * BE1 
1351:                 
1352:                CALL MTRXIN (APB, APBINV) 
1353:  
1354:                ARIBRJ =  LAMDAC1 * MATMUL(AE1,RI) + (1.D0 - LAMDAC1) * MATMUL(BE1,RJ) 
1355:                XC     =  MATMUL(APBINV, ARIBRJ) 
1356:                XCMRI  = XC - RI 
1357:                XCMRJ  = XC - RJ 
1358:                DF1DR  = - 2.D0 * LAMDAC1 * MATMUL(AE1,XCMRI) 
1359:  
1360:                D1ABEZ = MATMUL(DPI1RM,AEZR1(J1,:,:)) 
1361:                D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D1ABEZ))) 
1362:  
1363:                D2ABEZ = MATMUL(DPI2RM,AEZR1(J1,:,:)) 
1364:                D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D2ABEZ))) 
1365:  
1366:                D3ABEZ = MATMUL(DPI3RM,AEZR1(J1,:,:)) 
1367:                D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D3ABEZ))) 
1368:  
1369:                DF1PI1 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D1ABE,XCMRI)) 
1370:                DF1PI2 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D2ABE,XCMRI)) 
1371:                DF1PI3 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D3ABE,XCMRI)) 
1372:  
1373:                D1ABEZ = MATMUL(DPJ1RM,AEZR1(J2,:,:)) 
1374:                D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D1ABEZ))) 
1375:  
1376:                D2ABEZ = MATMUL(DPJ2RM,AEZR1(J2,:,:)) 
1377:                D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D2ABEZ))) 
1378:  
1379:                D3ABEZ = MATMUL(DPJ3RM,AEZR1(J2,:,:)) 
1380:                D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D3ABEZ)))  
1381:                 
1382:                DF1PJ1 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D1ABE,XCMRJ)) 
1383:                DF1PJ2 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D2ABE,XCMRJ)) 
1384:                DF1PJ3 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D3ABE,XCMRJ)) 
1385:  
1386:                RHO1   = PYSIGNOT / (ABSRIJ - ABSRIJ*SRTFI1 + PYSIGNOT) 
1387:                RHO1SQ = RHO1*RHO1 
1388:                RHO16  = RHO1SQ*RHO1SQ*RHO1SQ 
1389:                RHO112 = RHO16 * RHO16 
1390:  
1391:                FCTR1  = 0.5D0*ABSRIJ*SRTFI1/(FCNT1*PYSIGNOT) 
1392:                DG1DR  = (1.D0-SRTFI1)*NR/PYSIGNOT + FCTR1*DF1DR 
1393:                DVDF1  = -2.D0*RHO112*RHO1*FCTR1 
1394:  
1395:                IF (RADIFT) THEN 
1396:  
1397:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE2, BE2, RIJ, LAMDAC2, FMIN) 
1398:  
1399:                   FCNT2   = - FMIN 
1400:                   SRTFI2  = 1.D0 / DSQRT(FCNT2) 
1401:                   APB     = LAMDAC2 * AE2 + (1.D0 - LAMDAC2) * BE2 
1402:  
1403:                   CALL MTRXIN (APB, APBINV) 
1404:  
1405:                   ARIBRJ =  LAMDAC2 * MATMUL(AE2,RI) + (1.D0 - LAMDAC2) * MATMUL(BE2,RJ) 
1406:                   XC     =  MATMUL(APBINV, ARIBRJ) 
1407:                   XCMRI  = XC - RI 
1408:                   XCMRJ  = XC - RJ 
1409:                   DF2DR  = - 2.D0 * LAMDAC2 * MATMUL(AE2,XCMRI) 
1410:  
1411:                   RHO2   = PYSIGNOT / (ABSRIJ - ABSRIJ*SRTFI2 + PYSIGNOT) 
1412:                   RHO2SQ = RHO2*RHO2 
1413:                   RHO26  = RHO2SQ*RHO2SQ*RHO2SQ 
1414:                 
1415:                   FCTR2  = 0.5D0*ABSRIJ*SRTFI2/(FCNT2*PYSIGNOT) 
1416:                   DG2DR  = (1.D0-SRTFI2)*NR/PYSIGNOT+FCTR2*DF2DR 
1417:                   DVDF2  = RHO26*RHO2*FCTR2 
1418:  
1419:                   D1ABEZ = MATMUL(DPI1RM,AEZR2(J1,:,:)) 
1420:                   D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D1ABEZ))) 
1421:  
1422:                   D2ABEZ = MATMUL(DPI2RM,AEZR2(J1,:,:)) 
1423:                   D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D2ABEZ))) 
1424:  
1425:                   D3ABEZ = MATMUL(DPI3RM,AEZR2(J1,:,:)) 
1426:                   D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D3ABEZ))) 
1427:  
1428:                   DF2PI1 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D1ABE),XCMRI) 
1429:                   DF2PI2 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D2ABE),XCMRI) 
1430:                   DF2PI3 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D3ABE),XCMRI) 
1431:  
1432:                   D1ABEZ = MATMUL(DPJ1RM,AEZR2(J2,:,:)) 
1433:                   D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D1ABEZ))) 
1434:  
1435:                   D2ABEZ = MATMUL(DPJ2RM,AEZR2(J2,:,:)) 
1436:                   D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D2ABEZ))) 
1437:  
1438:                   D3ABEZ = MATMUL(DPJ3RM,AEZR2(J2,:,:)) 
1439:                   D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D3ABEZ))) 
1440:  
1441:                   DF2PJ1 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D1ABE),XCMRJ) 
1442:                   DF2PJ2 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D2ABE),XCMRJ) 
1443:                   DF2PJ3 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D3ABE),XCMRJ) 
1444:  
1445:                ELSE 
1446:  
1447:                   RHO2   = RHO1 
1448:                   RHO26  = RHO16 
1449:                   DG2DR  = DG1DR 
1450:                   DVDF2  = RHO26*RHO2*FCTR1 
1451:                   DF2PI1 = DF1PI1 
1452:                   DF2PI2 = DF1PI2 
1453:                   DF2PI3 = DF1PI3 
1454:                   DF2PJ1 = DF1PJ1 
1455:                   DF2PJ2 = DF1PJ2 
1456:                   DF2PJ3 = DF1PJ3 
1457:  
1458:                   SRTFI2 = SRTFI1 
1459:                   DF2DR  = DF1DR 
1460:                   RHO2SQ = RHO1SQ 
1461:  
1462:                ENDIF 
1463:  
1464:  
1465: ! Correction terms to the potential if we require a cutoff at rc.  
1466:                  SRTFI(1)=SRTFI1 
1467:                  SRTFI(2)=SRTFI2 
1468:                  DFP(1,1)=DF1PI1 
1469:                  DFP(1,2)=DF1PI2 
1470:                  DFP(1,3)=DF1PI3 
1471:                  DFP(1,4)=DF1PJ1 
1472:                  DFP(1,5)=DF1PJ2 
1473:                  DFP(1,6)=DF1PJ3 
1474:                  DFP(2,1)=DF2PI1 
1475:                  DFP(2,2)=DF2PI2 
1476:                  DFP(2,3)=DF2PI3 
1477:                  DFP(2,4)=DF2PJ1 
1478:                  DFP(2,5)=DF2PJ2 
1479:                  DFP(2,6)=DF2PJ3 
1480:                  DFDR(1,:)=DF1DR(:) 
1481:                  DFDR(2,:)=DF2DR(:) 
1482:                  LJ1(1)=RHO1 
1483:                  LJ1(2)=RHO2 
1484:  !     CALCULATE PY POTENTIAL ENERGY 
1485: 123 CONTINUE  
1486:                 VDUMMY=0.0D0 
1487:                IF(LJSITE) THEN   
1488:                 do k=1,maxinteractions ! k=1 -- interaction between repulsive primary 'apex' sites 
1489:                          ! k=2 and k=3 -- interaction between secondary and primary 'apex' sites 
1490:                          ! k=4 -- interaction between secondary 'apex' sites (normal LJ interaction)  
1491:                         ! trying to modify code to allow for binary systems.  
1492:                         ! apex site heights will be defined in absolute units, 
1493:                         ! hence PYA1bin(J1,1) etc. will be removed from below 
1494:  
1495:                    IF(k==1) THEN 
1496:                         DUMMY1=PSCALEFAC1vec(J1)!*PYA1bin(J1,1) 
1497:                         DUMMY2=PSCALEFAC1vec(J2)!*PYA1bin(J2,1) 
1498:                    ELSE IF(k==2) THEN 
1499:                         DUMMY1=PSCALEFAC1vec(J1)!*PYA1bin(J1,1) 
1500:                         DUMMY2=-PSCALEFAC2vec(J2)!*PYA1bin(J2,1) 
1501:                    ELSE IF(k==3) THEN 
1502:                         DUMMY1=-PSCALEFAC2vec(J1)!*PYA1bin(J1,1) 
1503:                         DUMMY2=PSCALEFAC1vec(J2)!*PYA1bin(J2,1) 
1504:                    ELSE 
1505:                         DUMMY1=-PSCALEFAC2vec(J1)!*PYA1bin(J1,1) 
1506:                         DUMMY2=-PSCALEFAC2vec(J2)!*PYA1bin(J2,1) 
1507:                    END IF 
1508:                         ! first particle 
1509:                         xlj(k,1,:)=RI+DUMMY1*MATMUL(RMI,vecsbf)    ! vecsbf: (1,0,0) in the body frame of ellipsoid 
1510:  
1511:                         ! second particle 
1512:                         xlj(k,2,:)=RJ+DUMMY2*MATMUL(RMJ,vecsbf) 
1513:  
1514:                         ! separation between the LJ sites 
1515:                         rlj2(k)=(xlj(k,2,1)-xlj(k,1,1))**2+(xlj(k,2,2)-xlj(k,1,2))**2+(xlj(k,2,3)-xlj(k,1,3))**2 
1516:                         rlj(k)=sqrt(rlj2(k)) 
1517:                         rljvec(k,1)=xlj(k,2,1)-xlj(k,1,1) 
1518:                         rljvec(k,2)=xlj(k,2,2)-xlj(k,1,2) 
1519:                         rljvec(k,3)=xlj(k,2,3)-xlj(k,1,3) 
1520:  
1521:                         DUMMY=1.0D0/rlj(k) 
1522:                         rljunitvec(k,:)=rljvec(k,:)*DUMMY !/rlj(k) 
1523:  
1524:                         drlj(k,1)=DUMMY*(xlj(k,2,1)-xlj(k,1,1))         !drlj/dx1 
1525:                         drlj(k,2)=DUMMY*(xlj(k,2,2)-xlj(k,1,2))         !drlj/dy1 
1526:                         drlj(k,3)=DUMMY*(xlj(k,2,3)-xlj(k,1,3))         !drlj/dz1 
1527:                         drlj(k,4)=-drlj(k,1)                               !drlj/dx2 
1528:                         drlj(k,5)=-drlj(k,2)                               !drlj/dy2 
1529:                         drlj(k,6)=-drlj(k,3)                               !drlj/dz2 
1530:                         drlj(k,7) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI1RM,vecsbf)) !drlj/dpx1 
1531:                         drlj(k,8) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI2RM,vecsbf)) !drlj/dpy1 
1532:                         drlj(k,9) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI3RM,vecsbf)) !drlj/dpz1 
1533:                         drlj(k,10) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ1RM,vecsbf)) !drlj/dpx2 
1534:                         drlj(k,11) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ2RM,vecsbf)) !drlj/dpy2 
1535:                         drlj(k,12) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ3RM,vecsbf)) !drlj/dpz2 
1536:  
1537:               ! interaction between the extra LJ sites: 
1538:                         LLJ(1,k)=sigma1(k)*DUMMY !/rlj(k) 
1539:                         LLJ(2,k)=LLJ(1,k)**2 
1540:                         LLJ(3,k)=LLJ(2,k)*LLJ(1,k) 
1541:                         LLJ(4,k)=LLJ(2,k)**2 
1542:                         LLJ(5,k)=LLJ(4,k)*LLJ(1,k) 
1543:                         LLJ(6,k)=LLJ(4,k)*LLJ(2,k) 
1544:                         LLJ(7,k)=LLJ(6,k)*LLJ(1,k) 
1545:                         LLJ(11,k)=LLJ(5,k)*LLJ(6,k) 
1546:                         LLJ(12,k)=LLJ(6,k)*LLJ(6,k) 
1547:  
1548:                             DO j=1,12 
1549:                                 dLLJ1(k,j) =-sigma1(k)*DUMMY*DUMMY*drlj(k,j) 
1550:                             END DO 
1551:  
1552:                  VDUMMY=VDUMMY+4.0D0*epsilon1(k,J1,J2)*term2(k)*(LLJ(12,k) - attr(k)*LLJ(6,k)) ! extra LJ sites (12-6) 
1553:               end do ! k=1,4 
1554:              END IF ! IF(LJSITE) 
1555:                 VDUMMY = VDUMMY + 4.0D0 * PYEPSNOT * (RHO112*PYADDREP(MJ2,MJ1) - RHO26*PYADDATT(MJ2,MJ1)) 
1556:  
1557:                ENERGY = ENERGY + VDUMMY 
1558:                VT(J1) = VT(J1) + VDUMMY 
1559:                VT(J2) = VT(J2) + VDUMMY        ! pair potentials 
1560:  
1561:             dVDUMMY(:) = 0.0D0 
1562:  
1563: !     CALCULATE GRADIENT 
1564:  
1565:              FIJ = 0.0D0 
1566:              TIJ = 0.0D0 
1567:              TJI = 0.0D0 
1568:              dVDUMMY(:)=0.0D0 
1569:  
1570:            IF(LJSITE) THEN 
1571:             do k=1,maxinteractions 
1572:              do j=1,3 
1573:                dVDUMMY(j) = dVDUMMY(j) + 4.0D0*epsilon1(k,J1,J2)*(12.0D0*LLJ(11,k)*dLLJ1(k,j))*term2(k) + & 
1574:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(12,k)*term3(k)*drlj(k,j)! extra LJ site derivatives, currently 
1575:                                                                                 ! only repulsive (LJ12) 
1576:                dVDUMMY(j) = dVDUMMY(j) - attr(k)*(4.0D0*epsilon1(k,J1,J2)*(6.0D0*LLJ(5,k)*dLLJ1(k,j))*term2(k) + & 
1577:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(6,k)*term3(k)*drlj(k,j)) ! attractive secondary apex site 
1578:                FIJ(j) = dVDUMMY(j) 
1579:              end do 
1580:              do j=7,12 
1581:                dVDUMMY(j) = dVDUMMY(j) + 4.0D0*epsilon1(k,J1,J2)*(12.0D0*LLJ(11,k)*dLLJ1(k,j))*term2(k) + & 
1582:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(12,k)*term3(k)*drlj(k,j)! extra LJ site derivatives, currently 
1583:                                                                                 ! only repulsive 
1584:                dVDUMMY(j) = dVDUMMY(j) - attr(k)*(4.0D0*epsilon1(k,J1,J2)*(6.0D0*LLJ(5,k)*dLLJ1(k,j))*term2(k) + & 
1585:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(6,k)*term3(k)*drlj(k,j)) ! attractive secondary apex site 
1586:              end do 
1587:             end do 
1588:            END IF !if LJSITE  
1589:                    FIJ    = FIJ + 24.0D0*PYEPSNOT*(2.D0*RHO112*RHO1*DG1DR*PYADDREP(MJ2,MJ1) - 1.0D0 * RHO26*RHO2*DG2DR*PYADDATT(MJ2,MJ1)) 
1590:                    TIJ(1) = DVDF1*DF1PI1*PYADDREP(MJ2,MJ1) + 1.0D0* DVDF2*DF2PI1*PYADDATT(MJ2,MJ1) 
1591:                    TIJ(2) = DVDF1*DF1PI2*PYADDREP(MJ2,MJ1) + 1.0D0*  DVDF2*DF2PI2*PYADDATT(MJ2,MJ1) 
1592:                    TIJ(3) = DVDF1*DF1PI3*PYADDREP(MJ2,MJ1) + 1.0D0*  DVDF2*DF2PI3*PYADDATT(MJ2,MJ1) 
1593:                    TJI(1) = DVDF1*DF1PJ1*PYADDREP(MJ2,MJ1) + 1.0D0*  DVDF2*DF2PJ1*PYADDATT(MJ2,MJ1) 
1594:                    TJI(2) = DVDF1*DF1PJ2*PYADDREP(MJ2,MJ1) + 1.0D0*  DVDF2*DF2PJ2*PYADDATT(MJ2,MJ1) 
1595:                    TJI(3) = DVDF1*DF1PJ3*PYADDREP(MJ2,MJ1) + 1.0D0*  DVDF2*DF2PJ3*PYADDATT(MJ2,MJ1) 
1596:                    TIJ = dVDUMMY(7:9) + 24.0D0 * PYEPSNOT * TIJ 
1597:                    TJI = dVDUMMY(10:12) + 24.0D0 * PYEPSNOT * TJI 
1598:  
1599:                G(J3-2:J3) = G(J3-2:J3) - FIJ 
1600:                G(J4-2:J4) = G(J4-2:J4) + FIJ 
1601:                G(J5-2:J5) = G(J5-2:J5) + TIJ 
1602:                G(J6-2:J6) = G(J6-2:J6) + TJI 
1603:  
1604: !     END INNER LOOP OVER PARTICLES 
1605: 124 CONTINUE 
1606:             ENDDO 
1607:  
1608: !     END OUTER LOOP OVER PARTICLES 
1609:  
1610:          ENDDO 
1611:         DO J1=1,NATOMS 
1612:                 J2=3*J1 
1613:                 IF(FROZEN(J1)) THEN 
1614:                         G(J2-2)=0.0D0 
1615:                         G(J2-1)=0.0D0 
1616:                         G(J2)=0.0D0 
1617:                 END IF 
1618:         END DO 
1619:         RETURN 
1620:       END SUBROUTINE PYGPERIODICADD2 
1621:  
1622:  
1623:  
1624: SUBROUTINE ECF(OVERLAPTEST,ECFvalue,x1,x2,y1,y2,z1,z2,px11,px12,py11,py12,pz11,pz12,a,b,c)1216: SUBROUTINE ECF(OVERLAPTEST,ECFvalue,x1,x2,y1,y2,z1,z2,px11,px12,py11,py12,pz11,pz12,a,b,c)
1625: 1217: 
1626: 1218: 
1627: !Calculates the value of the ellipsoid contact function1219: !Calculates the value of the ellipsoid contact function
1628: 1220: 
1629: !use commons1221: !use commons
1630: IMPLICIT NONE1222: IMPLICIT NONE
1631: 1223: 
1632: DOUBLE PRECISION, intent(in)       :: x1,x2,y1,y2,z1,z2,px11,px12,py11,py12,pz11,pz12,a,b,c1224: DOUBLE PRECISION, intent(in)       :: x1,x2,y1,y2,z1,z2,px11,px12,py11,py12,pz11,pz12,a,b,c
1633: DOUBLE PRECISION, intent(out)      :: ECFvalue1225: DOUBLE PRECISION, intent(out)      :: ECFvalue


r31913/io1.f 2017-02-15 08:30:16.924724926 +0000 r31912/io1.f 2017-02-15 08:30:18.092740538 +0000
630: 630: 
631:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' polarisableKZ water molecules '631:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' polarisableKZ water molecules '
632: 632: 
633:       ELSE IF (GAYBERNET) THEN633:       ELSE IF (GAYBERNET) THEN
634: 634: 
635:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' Gay-Berne atoms '635:          WRITE(MYUNIT,'(I4,A)') NATOMS/2,' Gay-Berne atoms '
636: 636: 
637:       ELSE IF (PYGPERIODICT) THEN637:       ELSE IF (PYGPERIODICT) THEN
638: 638: 
639: 639: 
640:          IF  (PYADDT.OR.PYADD2T) THEN640:          IF  (PYADDT) THEN
641:             IF (SORTT) THEN641:             IF (SORTT) THEN
642:                WRITE(MYUNIT,'(A)') 'Turning off SORT option for PYADD'642:                WRITE(MYUNIT,'(A)') 'Turning off SORT option for PYADD'
643:                SORTT=.FALSE.643:                SORTT=.FALSE.
644:             ENDIF644:             ENDIF
645:             WRITE(MYUNIT,'(I4,A)') NATOMS/2,' addressable ellipsoids'645:             WRITE(MYUNIT,'(I4,A)') NATOMS/2,' addressable ellipsoids'
646:          ELSE646:          ELSE
647:             WRITE(MYUNIT,'(I4,A)') NATOMS/2,' ellipsoids '647:             WRITE(MYUNIT,'(I4,A)') NATOMS/2,' ellipsoids '
648:          ENDIF648:          ENDIF
649: 649: 
650:       ELSE IF (PYT) THEN650:       ELSE IF (PYT) THEN


r31913/keywords.f 2017-02-15 08:30:17.152727974 +0000 r31912/keywords.f 2017-02-15 08:30:18.320743586 +0000
1167:       MLQT=.FALSE.1167:       MLQT=.FALSE.
1168:       MLQPROB=.FALSE.1168:       MLQPROB=.FALSE.
1169:       MLQDONE=.FALSE.1169:       MLQDONE=.FALSE.
1170:       MLQNORM=.FALSE.1170:       MLQNORM=.FALSE.
1171:       MLQLAMBDA=0.0D01171:       MLQLAMBDA=0.0D0
1172:       MLQSTART=11172:       MLQSTART=1
1173: 1173: 
1174:       LJADDT=.FALSE.1174:       LJADDT=.FALSE.
1175:       LJADD2T=.FALSE.1175:       LJADD2T=.FALSE.
1176:       PYADDT=.FALSE.1176:       PYADDT=.FALSE.
1177:       PYADD2T=.FALSE. 
1178: 1177: 
1179:       DUMPMQT=.FALSE.1178:       DUMPMQT=.FALSE.
1180: 1179: 
1181: 1180: 
1182: ! OPEP stuff1181: ! OPEP stuff
1183:       OPEPT = .FALSE.1182:       OPEPT = .FALSE.
1184:       OPEP_RNAT = .FALSE.1183:       OPEP_RNAT = .FALSE.
1185:              1184:              
1186:       CALL FILE_OPEN('data', DATA_UNIT, .FALSE.)1185:       CALL FILE_OPEN('data', DATA_UNIT, .FALSE.)
1187:       1186:       
6800:          LUNIT=GETUNIT()6799:          LUNIT=GETUNIT()
6801:          OPEN(LUNIT,FILE='epsilon',STATUS='OLD')6800:          OPEN(LUNIT,FILE='epsilon',STATUS='OLD')
6802:          IF (.NOT.ALLOCATED(PYADDEPS)) ALLOCATE(PYADDEPS(NADDTARGET,NADDTARGET))6801:          IF (.NOT.ALLOCATED(PYADDEPS)) ALLOCATE(PYADDEPS(NADDTARGET,NADDTARGET))
6803:          DO J1=1,NADDTARGET6802:          DO J1=1,NADDTARGET
6804:             DO J2=1,NADDTARGET6803:             DO J2=1,NADDTARGET
6805:                READ(LUNIT,*) PYADDEPS(J2,J1)6804:                READ(LUNIT,*) PYADDEPS(J2,J1)
6806:                WRITE(MYUNIT,'(2I6,G20.10)') J1,J2,PYADDEPS(J2,J1)6805:                WRITE(MYUNIT,'(2I6,G20.10)') J1,J2,PYADDEPS(J2,J1)
6807:             ENDDO6806:             ENDDO
6808:          ENDDO6807:          ENDDO
6809:          CLOSE(LUNIT)6808:          CLOSE(LUNIT)
6810:       ELSE IF (WORD.EQ.'PYADD2') THEN 
6811:          PYADD2T=.TRUE. 
6812:          CALL READI(NADDTARGET) 
6813:          WRITE(MYUNIT,'(A,I6)') 'keyword> Target cluster size is ',NADDTARGET 
6814:          IF (MOD(NATOMSALLOC/2,NADDTARGET).NE.0) THEN 
6815:             WRITE(MYUNIT,'(A,I6)') 'keyword> ERROR, target cluster size is not a factor of the number of PY particles ',  
6816:      &          NATOMSALLOC/2   
6817:             STOP 
6818:          ENDIF 
6819:          LUNIT=GETUNIT() 
6820:          OPEN(LUNIT,FILE='epsilon',STATUS='OLD') 
6821:          IF (.NOT.ALLOCATED(PYADDREP)) ALLOCATE(PYADDREP(NADDTARGET,NADDTARGET)) 
6822:          IF (.NOT.ALLOCATED(PYADDATT)) ALLOCATE(PYADDATT(NADDTARGET,NADDTARGET)) 
6823:          DO J1=1,NADDTARGET 
6824:             DO J2=1,NADDTARGET 
6825:                READ(LUNIT,*) PYADDREP(J2,J1), PYADDATT(J2,J1) 
6826:                WRITE(MYUNIT,'(2I6,2G20.10)') J1,J2,PYADDREP(J2,J1),PYADDATT(J2,J1) 
6827:             ENDDO 
6828:          ENDDO 
6829:          CLOSE(LUNIT) 
6830: 6809: 
6831:       ELSE IF (WORD .EQ.'PYGPERIODIC') THEN6810:       ELSE IF (WORD .EQ.'PYGPERIODIC') THEN
6832: 6811: 
6833:          PYGPERIODICT = .TRUE.6812:          PYGPERIODICT = .TRUE.
6834:          ELLIPSOIDT = .TRUE.6813:          ELLIPSOIDT = .TRUE.
6835:          RIGID = .TRUE.6814:          RIGID = .TRUE.
6836:          CALL READF(PYA1(1))6815:          CALL READF(PYA1(1))
6837:          CALL READF(PYA1(2))6816:          CALL READF(PYA1(2))
6838:          CALL READF(PYA1(3))6817:          CALL READF(PYA1(3))
6839:          CALL READF(PYA2(1))6818:          CALL READF(PYA2(1))


r31913/potential.f90 2017-02-15 08:30:17.424731611 +0000 r31912/potential.f90 2017-02-15 08:30:18.544746581 +0000
1126: 1126: 
1127:    ELSE IF (WATERKZT) THEN1127:    ELSE IF (WATERKZT) THEN
1128:       CALL WATERPKZ (X, GRAD, EREAL, GRADT)1128:       CALL WATERPKZ (X, GRAD, EREAL, GRADT)
1129: 1129: 
1130:    ELSE IF (GAYBERNET) THEN1130:    ELSE IF (GAYBERNET) THEN
1131:       CALL GAYBERNE(X, GRAD, EREAL, GRADT,.FALSE.)1131:       CALL GAYBERNE(X, GRAD, EREAL, GRADT,.FALSE.)
1132: 1132: 
1133:    ELSE IF (PYGPERIODICT .OR. PYBINARYT) THEN1133:    ELSE IF (PYGPERIODICT .OR. PYBINARYT) THEN
1134:       IF (PYADDT) THEN1134:       IF (PYADDT) THEN
1135:          CALL PYGPERIODICADD(X, GRAD, EREAL, GRADT)1135:          CALL PYGPERIODICADD(X, GRAD, EREAL, GRADT)
1136:       ELSEIF (PYADD2T) THEN 
1137:          CALL PYGPERIODICADD2(X, GRAD, EREAL, GRADT) 
1138:       ELSE1136:       ELSE
1139:          CALL PYGPERIODIC(X, GRAD, EREAL, GRADT)1137:          CALL PYGPERIODIC(X, GRAD, EREAL, GRADT)
1140:       ENDIF1138:       ENDIF
1141: 1139: 
1142:    ELSE IF (LJCAPSIDT) THEN1140:    ELSE IF (LJCAPSIDT) THEN
1143:       CALL LJCAPSIDMODEL(X, GRAD, EREAL, GRADT)1141:       CALL LJCAPSIDMODEL(X, GRAD, EREAL, GRADT)
1144: 1142: 
1145:    ELSE IF (STICKYT) THEN1143:    ELSE IF (STICKYT) THEN
1146:       CALL RADR(X, GRAD, EREAL, GRADT)1144:       CALL RADR(X, GRAD, EREAL, GRADT)
1147:       CALL STICKY(X, GRAD, EREAL, GRADT,.FALSE.)1145:       CALL STICKY(X, GRAD, EREAL, GRADT,.FALSE.)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0