hdiff output

r31914/gay-berne.f90 2017-02-15 10:30:13.756846478 +0000 r31913/gay-berne.f90 2017-02-15 10:30:14.648858378 +0000
1183:             end do1183:             end do
1184:            END IF1184:            END IF
1185:               IF (MJ1.EQ.MJ2) THEN ! same address - take just repulsion?1185:               IF (MJ1.EQ.MJ2) THEN ! same address - take just repulsion?
1186:                    FIJ    = FIJ + 24.0D0*2.D0*RHO112*RHO1*DG1DR*PYADDEPS(MJ2,MJ1)1186:                    FIJ    = FIJ + 24.0D0*2.D0*RHO112*RHO1*DG1DR*PYADDEPS(MJ2,MJ1)
1187:                    TIJ(1) = DVDF1*DF1PI11187:                    TIJ(1) = DVDF1*DF1PI1
1188:                    TIJ(2) = DVDF1*DF1PI21188:                    TIJ(2) = DVDF1*DF1PI2
1189:                    TIJ(3) = DVDF1*DF1PI31189:                    TIJ(3) = DVDF1*DF1PI3
1190:                    TJI(1) = DVDF1*DF1PJ11190:                    TJI(1) = DVDF1*DF1PJ1
1191:                    TJI(2) = DVDF1*DF1PJ21191:                    TJI(2) = DVDF1*DF1PJ2
1192:                    TJI(3) = DVDF1*DF1PJ31192:                    TJI(3) = DVDF1*DF1PJ3
1193:                    TIJ = dVDUMMY(7:9) + 24.0D0 * TIJ*PYADDEPS(MJ2,MJ1)1193:                    TIJ = dVDUMMY(7:9) + 24.0D0 * PYEPSNOT * TIJ*PYADDEPS(MJ2,MJ1)
1194:                    TJI = dVDUMMY(10:12) + 24.0D0 * TJI*PYADDEPS(MJ2,MJ1)1194:                    TJI = dVDUMMY(10:12) + 24.0D0 * PYEPSNOT * TJI*PYADDEPS(MJ2,MJ1)
1195:                ELSE1195:                ELSE
1196:                    FIJ    = FIJ + 24.0D0*(2.D0*RHO112*RHO1*DG1DR - RHO26*RHO2*DG2DR)*PYADDEPS(MJ2,MJ1)1196:                    FIJ    = FIJ + 24.0D0*(2.D0*RHO112*RHO1*DG1DR - RHO26*RHO2*DG2DR)*PYADDEPS(MJ2,MJ1)
1197:                    TIJ(1) = DVDF1*DF1PI1 + 1.0D0* DVDF2*DF2PI11197:                    TIJ(1) = DVDF1*DF1PI1 + 1.0D0* DVDF2*DF2PI1
1198:                    TIJ(2) = DVDF1*DF1PI2 + 1.0D0*  DVDF2*DF2PI21198:                    TIJ(2) = DVDF1*DF1PI2 + 1.0D0*  DVDF2*DF2PI2
1199:                    TIJ(3) = DVDF1*DF1PI3 + 1.0D0*  DVDF2*DF2PI31199:                    TIJ(3) = DVDF1*DF1PI3 + 1.0D0*  DVDF2*DF2PI3
1200:                    TJI(1) = DVDF1*DF1PJ1 + 1.0D0*  DVDF2*DF2PJ11200:                    TJI(1) = DVDF1*DF1PJ1 + 1.0D0*  DVDF2*DF2PJ1
1201:                    TJI(2) = DVDF1*DF1PJ2 + 1.0D0*  DVDF2*DF2PJ21201:                    TJI(2) = DVDF1*DF1PJ2 + 1.0D0*  DVDF2*DF2PJ2
1202:                    TJI(3) = DVDF1*DF1PJ3 + 1.0D0*  DVDF2*DF2PJ31202:                    TJI(3) = DVDF1*DF1PJ3 + 1.0D0*  DVDF2*DF2PJ3
1203:                    TIJ = dVDUMMY(7:9) + 24.0D0 * TIJ*PYADDEPS(MJ2,MJ1)1203:                    TIJ = dVDUMMY(7:9) + 24.0D0 * PYEPSNOT * TIJ*PYADDEPS(MJ2,MJ1)
1204:                    TJI = dVDUMMY(10:12) + 24.0D0 * TJI*PYADDEPS(MJ2,MJ1)1204:                    TJI = dVDUMMY(10:12) + 24.0D0 * PYEPSNOT * TJI*PYADDEPS(MJ2,MJ1)
1205:                ENDIF1205:                ENDIF
1206: 1206: 
1207:                G(J3-2:J3) = G(J3-2:J3) - FIJ1207:                G(J3-2:J3) = G(J3-2:J3) - FIJ
1208:                G(J4-2:J4) = G(J4-2:J4) + FIJ1208:                G(J4-2:J4) = G(J4-2:J4) + FIJ
1209:                G(J5-2:J5) = G(J5-2:J5) + TIJ1209:                G(J5-2:J5) = G(J5-2:J5) + TIJ
1210:                G(J6-2:J6) = G(J6-2:J6) + TJI1210:                G(J6-2:J6) = G(J6-2:J6) + TJI
1211: 1211: 
1212: !     END INNER LOOP OVER PARTICLES1212: !     END INNER LOOP OVER PARTICLES
1213: 124 CONTINUE1213: 124 CONTINUE
1214:             ENDDO1214:             ENDDO
1216: !     END OUTER LOOP OVER PARTICLES1216: !     END OUTER LOOP OVER PARTICLES
1217: 1217: 
1218:          ENDDO1218:          ENDDO
1219:          ENERGY = PYEPSNOT*ENERGY1219:          ENERGY = PYEPSNOT*ENERGY
1220:          G      = PYEPSNOT*G1220:          G      = PYEPSNOT*G
1221: 1221: 
1222:         RETURN1222:         RETURN
1223:       END SUBROUTINE PYGPERIODICADD1223:       END SUBROUTINE PYGPERIODICADD
1224: 1224: 
1225: 1225: 
1226: ! PY potential, dc430's implementation 
1227: ! with PBC and continuous cutoff added 
1228: ! plus extra LJ site 
1229: ! 
1230: ! This version is addressable with interactions scaled by the input epsilon matrix PYADDEPS  
1231: ! with separate rep and att entries 
1232: ! 
1233:       SUBROUTINE PYGPERIODICADD2 (X, G, ENERGY, GTEST, STEST) 
1234:  
1235:        use key, only: PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PCUTOFF,PARAMONOVCUTOFF,& 
1236:                 &       pya1bin,pya2bin,pysignot,pyepsnot,radift,LJSITE,BLJSITE,PEPSILON1,& 
1237:                 &       PSCALEFAC1,PSCALEFAC2,MAXINTERACTIONS,PYBINARYT,PYBINARYTYPE1,VT,PEPSILONATTR,PSIGMAATTR,PYADDREP, PYADDATT, NADDTARGET 
1238:        use commons, only: natoms 
1239:        use pymodule 
1240:       IMPLICIT NONE 
1241:  
1242:       DOUBLE PRECISION :: X(3*NATOMS), G(3*NATOMS) 
1243:       DOUBLE PRECISION :: RI(3), RJ(3), RIJ(3), NR(3), RIJSQ, ABSRIJ, P(3), THETA, THETA2, CT, ST 
1244:       DOUBLE PRECISION :: AE1(3,3), BE1(3,3), AE2(3,3), BE2(3,3), RM(3,3), APB(3,3), APBINV(3,3) 
1245:       DOUBLE PRECISION :: FCNT1, FCNT2, SRTFI1, SRTFI2, SRTFI(2), FMIN, LAMDAC1, LAMDAC2, ENERGY 
1246:       DOUBLE PRECISION :: RHO1, RHO1SQ, RHO16, RHO112, RHO2, RHO2SQ, RHO26 
1247:       DOUBLE PRECISION :: FCTR1, FCTR2, DVDF1, DVDF2  
1248:       DOUBLE PRECISION :: DF1PI1, DF1PI2, DF1PI3, DF2PI1, DF2PI2, DF2PI3 
1249:       DOUBLE PRECISION :: DF1PJ1, DF1PJ2, DF1PJ3, DF2PJ1, DF2PJ2, DF2PJ3  
1250:       DOUBLE PRECISION :: RMI(3,3), RMJ(3,3), E(3,3), ESQ(3,3), DE1(3,3), DE2(3,3), DE3(3,3) 
1251:       DOUBLE PRECISION :: DPI1RM(3,3), DPI2RM(3,3), DPI3RM(3,3), DPJ1RM(3,3), DPJ2RM(3,3), DPJ3RM(3,3) 
1252:       DOUBLE PRECISION :: DF1DR(3), DF2DR(3), DG1DR(3), DG2DR(3) 
1253:       DOUBLE PRECISION :: ARIBRJ(3), XC(3), XCMRI(3), XCMRJ(3), FIJ(3), TIJ(3), TJI(3) 
1254:       double precision :: D1ABEZ(3,3), D2ABEZ(3,3), D3ABEZ(3,3), D1ABE(3,3), D2ABE(3,3), D3ABE(3,3)  
1255:       LOGICAL          :: GTEST 
1256:  
1257:  
1258: ! sf344 additions 
1259:  
1260:       DOUBLE PRECISION :: CLJ1(2),CLJ2(2),CLJ3(2),CLJ4(2),CLJ5(2),CLJ6(2),CLJ7(2),CLJ8(2),CLJ11(2),CLJ12(2), & 
1261:                           & CLJ13(2),CLJ14(2),DFDR(2,3),DFP(2,6),dr(6),dCLJ1(2,12),dVDUMMY(12),LJ1(2),DUMMY,DUMMY1,DUMMY2,& 
1262:                           & dLJ1(2,12),VDUMMY 
1263:       DOUBLE PRECISION :: term2(maxinteractions),term3(maxinteractions), & 
1264:                           & xlj(maxinteractions,2,3),rljvec(maxinteractions,3),rljunitvec(maxinteractions,3),& 
1265:                           & drlj(maxinteractions,12),rlj(maxinteractions),rlj2(maxinteractions) 
1266:       DOUBLE PRECISION :: LLJ(12,maxinteractions), dLLJ1(maxinteractions,12) 
1267:       INTEGER          :: k, MJ1, MJ2 
1268:       LOGICAL          :: STEST 
1269:       DOUBLE PRECISION :: cosangle, sinangle, IMAT(3,3), tildematrix(3,3), rotmat(3,3), crossvector(3),tempcrd(6) 
1270:      VT(1:NATOMS/2)=0.0D0 
1271:      term2(:)=1.0D0 
1272:      term3(:)=0.0D0 
1273:  
1274:          ENERGY = 0.D0 
1275:          G(:)   = 0.D0 
1276:          
1277:         DO J1=1, REALNATOMS 
1278:  
1279:             J3      = 3*J1 
1280:             J5      = OFFSET + J3 
1281:             RI      = X(J3-2:J3) 
1282:             P       = X(J5-2:J5) 
1283:             angle=sqrt(dot_product(P,P)) 
1284:             if(angle>twopi) then 
1285: ! normalise angle-axis coordinates 
1286:                 X(J5-2:J5)=X(J5-2:J5)/angle 
1287:                 do 
1288:                   angle=angle-twopi 
1289:                   if(angle<2*pi) exit 
1290:                 end do 
1291: ! multiply with new angle 
1292:                 X(J5-2:J5)=X(J5-2:J5)*angle 
1293:             end if 
1294:  
1295:             CALL RMDRVT(P, RMIvec(J1,:,:), DPI1RMvec(J1,:,:), DPI2RMvec(J1,:,:), DPI3RMvec(J1,:,:), .TRUE.) 
1296:  
1297:         END DO 
1298:  
1299:          DO J1 = 1, REALNATOMS - 1 
1300:  
1301:             MJ1=MOD(J1-1,NADDTARGET)+1 
1302:             J3      = 3*J1 
1303:             J5      = OFFSET + J3 
1304:             RI      = X(J3-2:J3) 
1305:             P       = X(J5-2:J5) 
1306: !     ROTATION MATRIX 
1307:  
1308: !            CALL RMDRVT(P, RMI, DPI1RM, DPI2RM, DPI3RM, GTEST) 
1309:             RMI(:,:)=RMIvec(J1,:,:) 
1310:             DPI1RM(:,:)=DPI1RMvec(J1,:,:) 
1311:             DPI2RM(:,:)=DPI2RMvec(J1,:,:) 
1312:             DPI3RM(:,:)=DPI3RMvec(J1,:,:) 
1313:  
1314:             AE1 = MATMUL(RMI,(MATMUL(AEZR1(J1,:,:),(TRANSPOSE(RMI))))) 
1315:  
1316:             IF (RADIFT) THEN 
1317:  
1318:                AE2 = MATMUL(RMI,(MATMUL(AEZR2(J1,:,:),(TRANSPOSE(RMI)))))          
1319:  
1320:             ENDIF 
1321:  
1322: !     BEGIN INNER LOOP OVER PARTICLES 
1323:  
1324:             DO J2 = J1 + 1, REALNATOMS 
1325:  
1326:                MJ2=MOD(J2-1,NADDTARGET)+1 
1327:                J4     = 3*J2 
1328:                J6     = OFFSET + J4 
1329:                RJ     = X(J4-2:J4)  
1330:                P      = X(J6-2:J6) 
1331:  
1332: !     ROTATION MATRIX 
1333:  
1334: !               CALL RMDRVT(P, RMJ, DPJ1RM, DPJ2RM, DPJ3RM, GTEST)                
1335:                RMJ(:,:)=RMIvec(J2,:,:) 
1336:                DPJ1RM(:,:)=DPI1RMvec(J2,:,:) 
1337:                DPJ2RM(:,:)=DPI2RMvec(J2,:,:) 
1338:                DPJ3RM(:,:)=DPI3RMvec(J2,:,:) 
1339:       
1340:                BE1 = MATMUL(RMJ,(MATMUL(AEZR1(J2,:,:),(TRANSPOSE(RMJ))))) 
1341:  
1342:                IF (RADIFT) THEN 
1343:     
1344:                   BE2 = MATMUL(RMJ,(MATMUL(AEZR2(J2,:,:),(TRANSPOSE(RMJ))))) 
1345:  
1346:                ENDIF 
1347:  
1348: !     CALCULATE SEPARATION 
1349:  
1350:                RIJ    = RI - RJ 
1351:                RIJSQ  = DOT_PRODUCT(RIJ,RIJ) 
1352:                ABSRIJ = DSQRT(RIJSQ) 
1353:                NR     = RIJ / ABSRIJ 
1354:  
1355:                 IF(PARAMONOVCUTOFF.AND.RIJSQ>cut) GOTO 124 
1356:  
1357: !     CALCULATE ECF 
1358:  
1359:                CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE1, BE1, RIJ, LAMDAC1, FMIN) 
1360:  
1361:                FCNT1   = - FMIN 
1362:                SRTFI1  = 1.D0 / DSQRT(FCNT1) 
1363:                APB     = LAMDAC1 * AE1 + (1.D0 - LAMDAC1) * BE1 
1364:                 
1365:                CALL MTRXIN (APB, APBINV) 
1366:  
1367:                ARIBRJ =  LAMDAC1 * MATMUL(AE1,RI) + (1.D0 - LAMDAC1) * MATMUL(BE1,RJ) 
1368:                XC     =  MATMUL(APBINV, ARIBRJ) 
1369:                XCMRI  = XC - RI 
1370:                XCMRJ  = XC - RJ 
1371:                DF1DR  = - 2.D0 * LAMDAC1 * MATMUL(AE1,XCMRI) 
1372:  
1373:                D1ABEZ = MATMUL(DPI1RM,AEZR1(J1,:,:)) 
1374:                D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D1ABEZ))) 
1375:  
1376:                D2ABEZ = MATMUL(DPI2RM,AEZR1(J1,:,:)) 
1377:                D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D2ABEZ))) 
1378:  
1379:                D3ABEZ = MATMUL(DPI3RM,AEZR1(J1,:,:)) 
1380:                D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D3ABEZ))) 
1381:  
1382:                DF1PI1 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D1ABE,XCMRI)) 
1383:                DF1PI2 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D2ABE,XCMRI)) 
1384:                DF1PI3 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D3ABE,XCMRI)) 
1385:  
1386:                D1ABEZ = MATMUL(DPJ1RM,AEZR1(J2,:,:)) 
1387:                D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D1ABEZ))) 
1388:  
1389:                D2ABEZ = MATMUL(DPJ2RM,AEZR1(J2,:,:)) 
1390:                D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D2ABEZ))) 
1391:  
1392:                D3ABEZ = MATMUL(DPJ3RM,AEZR1(J2,:,:)) 
1393:                D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D3ABEZ)))  
1394:                 
1395:                DF1PJ1 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D1ABE,XCMRJ)) 
1396:                DF1PJ2 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D2ABE,XCMRJ)) 
1397:                DF1PJ3 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D3ABE,XCMRJ)) 
1398:  
1399:                RHO1   = PYSIGNOT / (ABSRIJ - ABSRIJ*SRTFI1 + PYSIGNOT) 
1400:                RHO1SQ = RHO1*RHO1 
1401:                RHO16  = RHO1SQ*RHO1SQ*RHO1SQ 
1402:                RHO112 = RHO16 * RHO16 
1403:  
1404:                FCTR1  = 0.5D0*ABSRIJ*SRTFI1/(FCNT1*PYSIGNOT) 
1405:                DG1DR  = (1.D0-SRTFI1)*NR/PYSIGNOT + FCTR1*DF1DR 
1406:                DVDF1  = -2.D0*RHO112*RHO1*FCTR1 
1407:  
1408:                IF (RADIFT) THEN 
1409:  
1410:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE2, BE2, RIJ, LAMDAC2, FMIN) 
1411:  
1412:                   FCNT2   = - FMIN 
1413:                   SRTFI2  = 1.D0 / DSQRT(FCNT2) 
1414:                   APB     = LAMDAC2 * AE2 + (1.D0 - LAMDAC2) * BE2 
1415:  
1416:                   CALL MTRXIN (APB, APBINV) 
1417:  
1418:                   ARIBRJ =  LAMDAC2 * MATMUL(AE2,RI) + (1.D0 - LAMDAC2) * MATMUL(BE2,RJ) 
1419:                   XC     =  MATMUL(APBINV, ARIBRJ) 
1420:                   XCMRI  = XC - RI 
1421:                   XCMRJ  = XC - RJ 
1422:                   DF2DR  = - 2.D0 * LAMDAC2 * MATMUL(AE2,XCMRI) 
1423:  
1424:                   RHO2   = PYSIGNOT / (ABSRIJ - ABSRIJ*SRTFI2 + PYSIGNOT) 
1425:                   RHO2SQ = RHO2*RHO2 
1426:                   RHO26  = RHO2SQ*RHO2SQ*RHO2SQ 
1427:                 
1428:                   FCTR2  = 0.5D0*ABSRIJ*SRTFI2/(FCNT2*PYSIGNOT) 
1429:                   DG2DR  = (1.D0-SRTFI2)*NR/PYSIGNOT+FCTR2*DF2DR 
1430:                   DVDF2  = RHO26*RHO2*FCTR2 
1431:  
1432:                   D1ABEZ = MATMUL(DPI1RM,AEZR2(J1,:,:)) 
1433:                   D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D1ABEZ))) 
1434:  
1435:                   D2ABEZ = MATMUL(DPI2RM,AEZR2(J1,:,:)) 
1436:                   D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D2ABEZ))) 
1437:  
1438:                   D3ABEZ = MATMUL(DPI3RM,AEZR2(J1,:,:)) 
1439:                   D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D3ABEZ))) 
1440:  
1441:                   DF2PI1 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D1ABE),XCMRI) 
1442:                   DF2PI2 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D2ABE),XCMRI) 
1443:                   DF2PI3 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D3ABE),XCMRI) 
1444:  
1445:                   D1ABEZ = MATMUL(DPJ1RM,AEZR2(J2,:,:)) 
1446:                   D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D1ABEZ))) 
1447:  
1448:                   D2ABEZ = MATMUL(DPJ2RM,AEZR2(J2,:,:)) 
1449:                   D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D2ABEZ))) 
1450:  
1451:                   D3ABEZ = MATMUL(DPJ3RM,AEZR2(J2,:,:)) 
1452:                   D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D3ABEZ))) 
1453:  
1454:                   DF2PJ1 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D1ABE),XCMRJ) 
1455:                   DF2PJ2 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D2ABE),XCMRJ) 
1456:                   DF2PJ3 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D3ABE),XCMRJ) 
1457:  
1458:                ELSE 
1459:  
1460:                   RHO2   = RHO1 
1461:                   RHO26  = RHO16 
1462:                   DG2DR  = DG1DR 
1463:                   DVDF2  = RHO26*RHO2*FCTR1 
1464:                   DF2PI1 = DF1PI1 
1465:                   DF2PI2 = DF1PI2 
1466:                   DF2PI3 = DF1PI3 
1467:                   DF2PJ1 = DF1PJ1 
1468:                   DF2PJ2 = DF1PJ2 
1469:                   DF2PJ3 = DF1PJ3 
1470:  
1471:                   SRTFI2 = SRTFI1 
1472:                   DF2DR  = DF1DR 
1473:                   RHO2SQ = RHO1SQ 
1474:  
1475:                ENDIF 
1476:  
1477:  
1478: ! Correction terms to the potential if we require a cutoff at rc.  
1479:                  SRTFI(1)=SRTFI1 
1480:                  SRTFI(2)=SRTFI2 
1481:                  DFP(1,1)=DF1PI1 
1482:                  DFP(1,2)=DF1PI2 
1483:                  DFP(1,3)=DF1PI3 
1484:                  DFP(1,4)=DF1PJ1 
1485:                  DFP(1,5)=DF1PJ2 
1486:                  DFP(1,6)=DF1PJ3 
1487:                  DFP(2,1)=DF2PI1 
1488:                  DFP(2,2)=DF2PI2 
1489:                  DFP(2,3)=DF2PI3 
1490:                  DFP(2,4)=DF2PJ1 
1491:                  DFP(2,5)=DF2PJ2 
1492:                  DFP(2,6)=DF2PJ3 
1493:                  DFDR(1,:)=DF1DR(:) 
1494:                  DFDR(2,:)=DF2DR(:) 
1495:                  LJ1(1)=RHO1 
1496:                  LJ1(2)=RHO2 
1497:  !     CALCULATE PY POTENTIAL ENERGY 
1498: 123 CONTINUE  
1499:         VDUMMY = 0.0D0 
1500:                   
1501:                IF(LJSITE) THEN   
1502:                 do k=1,maxinteractions ! k=1 -- interaction between repulsive primary 'apex' sites 
1503:                          ! k=2 and k=3 -- interaction between secondary and primary 'apex' sites  
1504:                          ! k=4 -- interaction between secondary 'apex' sites (normal LJ interaction)  
1505:                         ! trying to modify code to allow for binary systems.  
1506:                         ! apex site heights will be defined in absolute units, 
1507:                         ! hence PYA1bin(J1,1) etc. will be removed from below 
1508:  
1509:                    IF(k==1) THEN 
1510:                         DUMMY1=PSCALEFAC1vec(J1)!*PYA1bin(J1,1) 
1511:                         DUMMY2=PSCALEFAC1vec(J2)!*PYA1bin(J2,1) 
1512:                    ELSE IF(k==2) THEN 
1513:                         DUMMY1=PSCALEFAC1vec(J1)!*PYA1bin(J1,1) 
1514:                         DUMMY2=-PSCALEFAC2vec(J2)!*PYA1bin(J2,1) 
1515:                    ELSE IF(k==3) THEN 
1516:                         DUMMY1=-PSCALEFAC2vec(J1)!*PYA1bin(J1,1) 
1517:                         DUMMY2=PSCALEFAC1vec(J2)!*PYA1bin(J2,1) 
1518:                    ELSE  
1519:                         DUMMY1=-PSCALEFAC2vec(J1)!*PYA1bin(J1,1) 
1520:                         DUMMY2=-PSCALEFAC2vec(J2)!*PYA1bin(J2,1) 
1521:                    END IF 
1522:                         ! first particle 
1523:                         xlj(k,1,:)=RI+DUMMY1*MATMUL(RMI,vecsbf)    ! vecsbf: (1,0,0) in the body frame of ellipsoid 
1524:  
1525:                         ! second particle 
1526:                         xlj(k,2,:)=RJ+DUMMY2*MATMUL(RMJ,vecsbf) 
1527:  
1528:                         ! separation between the LJ sites 
1529:                         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 
1530:                         rlj(k)=sqrt(rlj2(k)) 
1531:                         rljvec(k,1)=xlj(k,2,1)-xlj(k,1,1) 
1532:                         rljvec(k,2)=xlj(k,2,2)-xlj(k,1,2) 
1533:                         rljvec(k,3)=xlj(k,2,3)-xlj(k,1,3) 
1534:  
1535:                         DUMMY=1.0D0/rlj(k) 
1536:                         rljunitvec(k,:)=rljvec(k,:)*DUMMY !/rlj(k) 
1537:  
1538:                         drlj(k,1)=DUMMY*(xlj(k,2,1)-xlj(k,1,1))         !drlj/dx1 
1539:                         drlj(k,2)=DUMMY*(xlj(k,2,2)-xlj(k,1,2))         !drlj/dy1 
1540:                         drlj(k,3)=DUMMY*(xlj(k,2,3)-xlj(k,1,3))         !drlj/dz1 
1541:                         drlj(k,4)=-drlj(k,1)                               !drlj/dx2 
1542:                         drlj(k,5)=-drlj(k,2)                               !drlj/dy2 
1543:                         drlj(k,6)=-drlj(k,3)                               !drlj/dz2 
1544:                         drlj(k,7) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI1RM,vecsbf)) !drlj/dpx1 
1545:                         drlj(k,8) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI2RM,vecsbf)) !drlj/dpy1 
1546:                         drlj(k,9) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI3RM,vecsbf)) !drlj/dpz1 
1547:                         drlj(k,10) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ1RM,vecsbf)) !drlj/dpx2 
1548:                         drlj(k,11) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ2RM,vecsbf)) !drlj/dpy2 
1549:                         drlj(k,12) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ3RM,vecsbf)) !drlj/dpz2 
1550:  
1551:               ! interaction between the extra LJ sites: 
1552:                         LLJ(1,k)=sigma1(k)*DUMMY !/rlj(k) 
1553:                         LLJ(2,k)=LLJ(1,k)**2 
1554:                         LLJ(3,k)=LLJ(2,k)*LLJ(1,k) 
1555:                         LLJ(4,k)=LLJ(2,k)**2 
1556:                         LLJ(5,k)=LLJ(4,k)*LLJ(1,k) 
1557:                         LLJ(6,k)=LLJ(4,k)*LLJ(2,k) 
1558:                         LLJ(7,k)=LLJ(6,k)*LLJ(1,k) 
1559:                         LLJ(11,k)=LLJ(5,k)*LLJ(6,k) 
1560:                         LLJ(12,k)=LLJ(6,k)*LLJ(6,k) 
1561:  
1562:                             DO j=1,12 
1563:                                 dLLJ1(k,j) =-sigma1(k)*DUMMY*DUMMY*drlj(k,j) 
1564:                             END DO 
1565:  
1566:                  VDUMMY=VDUMMY+4.0D0*epsilon1(k,J1,J2)*term2(k)*(LLJ(12,k) - attr(k)*LLJ(6,k)) ! extra LJ sites 
1567:               end do ! k=1,4 
1568:              END IF ! IF(LJSITE) 
1569:                 VDUMMY = VDUMMY + 4.0D0 * (RHO112*PYADDREP(MJ2,MJ1) - RHO26*PYADDATT(MJ2,MJ1)) 
1570:  
1571:                ENERGY = ENERGY + VDUMMY 
1572:                VT(J1) = VT(J1) + VDUMMY 
1573:                VT(J2) = VT(J2) + VDUMMY        ! pair potentials 
1574:  
1575:             dVDUMMY(:) = 0.0D0 
1576:  
1577: !     CALCULATE GRADIENT 
1578:              FIJ = 0.0D0 
1579:              TIJ = 0.0D0 
1580:              TJI = 0.0D0 
1581:              dVDUMMY(:)=0.0D0 
1582:  
1583:            IF(LJSITE) THEN 
1584:             do k=1,maxinteractions 
1585:              do j=1,3 
1586:                dVDUMMY(j) = dVDUMMY(j) + 4.0D0*epsilon1(k,J1,J2)*(12.0D0*LLJ(11,k)*dLLJ1(k,j))*term2(k) + & 
1587:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(12,k)*term3(k)*drlj(k,j)! extra LJ site derivatives, currently only repulsive 
1588:                dVDUMMY(j) = dVDUMMY(j) - attr(k)*(4.0D0*epsilon1(k,J1,J2)*(6.0D0*LLJ(5,k)*dLLJ1(k,j))*term2(k) + & 
1589:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(6,k)*term3(k)*drlj(k,j)) ! attractive secondary apex site 
1590:                FIJ(j) = dVDUMMY(j) 
1591:              end do 
1592:              do j=7,12 
1593:                dVDUMMY(j) = dVDUMMY(j) + 4.0D0*epsilon1(k,J1,J2)*(12.0D0*LLJ(11,k)*dLLJ1(k,j))*term2(k) + & 
1594:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(12,k)*term3(k)*drlj(k,j)! extra LJ site derivatives, currently only repulsive 
1595:                dVDUMMY(j) = dVDUMMY(j) - attr(k)*(4.0D0*epsilon1(k,J1,J2)*(6.0D0*LLJ(5,k)*dLLJ1(k,j))*term2(k) + & 
1596:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(6,k)*term3(k)*drlj(k,j)) ! attractive secondary apex site 
1597:              end do 
1598:             end do 
1599:            END IF 
1600:  
1601:                    FIJ    = FIJ + 24.0D0*(2.D0*RHO112*RHO1*DG1DR*PYADDREP(MJ2,MJ1) - 1.0D0 * RHO26*RHO2*DG2DR*PYADDATT(MJ2,MJ1)) 
1602:                    TIJ(1) = DVDF1*DF1PI1*PYADDREP(MJ2,MJ1) + 1.0D0* DVDF2*DF2PI1*PYADDATT(MJ2,MJ1) 
1603:                    TIJ(2) = DVDF1*DF1PI2*PYADDREP(MJ2,MJ1) + 1.0D0*  DVDF2*DF2PI2*PYADDATT(MJ2,MJ1) 
1604:                    TIJ(3) = DVDF1*DF1PI3*PYADDREP(MJ2,MJ1) + 1.0D0*  DVDF2*DF2PI3*PYADDATT(MJ2,MJ1) 
1605:                    TJI(1) = DVDF1*DF1PJ1*PYADDREP(MJ2,MJ1) + 1.0D0*  DVDF2*DF2PJ1*PYADDATT(MJ2,MJ1) 
1606:                    TJI(2) = DVDF1*DF1PJ2*PYADDREP(MJ2,MJ1) + 1.0D0*  DVDF2*DF2PJ2*PYADDATT(MJ2,MJ1) 
1607:                    TJI(3) = DVDF1*DF1PJ3*PYADDREP(MJ2,MJ1) + 1.0D0*  DVDF2*DF2PJ3*PYADDATT(MJ2,MJ1) 
1608:                    TIJ = dVDUMMY(7:9) + 24.0D0 * TIJ 
1609:                    TJI = dVDUMMY(10:12) + 24.0D0 * TJI 
1610:  
1611:                G(J3-2:J3) = G(J3-2:J3) - FIJ 
1612:                G(J4-2:J4) = G(J4-2:J4) + FIJ 
1613:                G(J5-2:J5) = G(J5-2:J5) + TIJ 
1614:                G(J6-2:J6) = G(J6-2:J6) + TJI 
1615:  
1616: !     END INNER LOOP OVER PARTICLES 
1617: 124 CONTINUE 
1618:             ENDDO 
1619:  
1620: !     END OUTER LOOP OVER PARTICLES 
1621:  
1622:          ENDDO 
1623:          ENERGY = PYEPSNOT*ENERGY 
1624:          G      = PYEPSNOT*G 
1625:  
1626:         RETURN 
1627:       END SUBROUTINE PYGPERIODICADD2 
1628:  
1629:  
1630: SUBROUTINE PARAMONOVNUMFIRSTDER(OLDX,STEST)1226: SUBROUTINE PARAMONOVNUMFIRSTDER(OLDX,STEST)
1631: use commons, only : natoms1227: use commons, only : natoms
1632: !use key1228: !use key
1633: use modhess1229: use modhess
1634: implicit none1230: implicit none
1635: 1231: 
1636: DOUBLE PRECISION  :: V(3*NATOMS),ENERGY,X(3*NATOMS),NUMGRAD(3*NATOMS),TRUEGRAD(3*NATOMS),&1232: DOUBLE PRECISION  :: V(3*NATOMS),ENERGY,X(3*NATOMS),NUMGRAD(3*NATOMS),TRUEGRAD(3*NATOMS),&
1637:         &            OLDX(3*NATOMS),ETEMP(2,3*NATOMS),ksi1233:         &            OLDX(3*NATOMS),ETEMP(2,3*NATOMS),ksi
1638: INTEGER    :: i,j1234: INTEGER    :: i,j
1639: LOGICAL    :: GTEST,STEST1235: LOGICAL    :: GTEST,STEST
1856:          VTEMP(2,:)=V(:)1452:          VTEMP(2,:)=V(:)
1857: 1453: 
1858:                 DO j=i,3*NATOMS1454:                 DO j=i,3*NATOMS
1859:                         HESS(i,j)=(VTEMP(2,j)-VTEMP(1,j))/(2.0D0*ksi)1455:                         HESS(i,j)=(VTEMP(2,j)-VTEMP(1,j))/(2.0D0*ksi)
1860:                         HESS(j,i)=HESS(i,j)1456:                         HESS(j,i)=HESS(i,j)
1861:                 END DO1457:                 END DO
1862: END DO1458: END DO
1863: 1459: 
1864: END SUBROUTINE PYGPERIODICSECDERADD1460: END SUBROUTINE PYGPERIODICSECDERADD
1865: 1461: 
1866:  
1867: SUBROUTINE PYGPERIODICSECDERADD2(OLDX,STEST) 
1868: use commons 
1869: use modhess 
1870: implicit none 
1871:  
1872: DOUBLE PRECISION  :: V(3*NATOMS),EGB,X(3*NATOMS),OLDX(3*NATOMS),VTEMP(2,3*NATOMS),ksi 
1873: INTEGER    :: i,j 
1874: LOGICAL    :: GTEST,STEST 
1875:  
1876: ksi=0.00001D0 
1877: X(:)=OLDX(:) 
1878:  
1879: IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS,3*NATOMS)) 
1880: HESS(:,:)=0.0D0 
1881: V(:)=0.0D0 
1882: VTEMP(:,:)=0.0D0 
1883:  
1884: DO i=1,3*NATOMS 
1885:  
1886:          X(i)=X(i)-ksi 
1887:   
1888:          CALL PYGPERIODICADD2(X,V,EGB,.TRUE.,.FALSE.) 
1889:          VTEMP(1,:)=V(:) 
1890:   
1891:          X(i)=X(i)+2.0D0*ksi 
1892:  
1893:          CALL PYGPERIODICADD2(X,V,EGB,.TRUE.,.FALSE.) 
1894:  
1895:          VTEMP(2,:)=V(:) 
1896:  
1897:                 DO j=i,3*NATOMS 
1898:                         HESS(i,j)=(VTEMP(2,j)-VTEMP(1,j))/(2.0D0*ksi) 
1899:                         HESS(j,i)=HESS(i,j) 
1900:                 END DO 
1901: END DO 
1902:  
1903: END SUBROUTINE PYGPERIODICSECDERADD2 
1904:  
1905: SUBROUTINE EllipsoidsAAtoPolar(px1,py1,pz1,alpha,beta,gamma,alphadeg,betadeg,gammadeg)1462: SUBROUTINE EllipsoidsAAtoPolar(px1,py1,pz1,alpha,beta,gamma,alphadeg,betadeg,gammadeg)
1906: ! converts angle-axis coordinates px, py, pz to "polar-like" angles alpha and beta1463: ! converts angle-axis coordinates px, py, pz to "polar-like" angles alpha and beta
1907: ! px=cos(alpha)*cos(alpha)1464: ! px=cos(alpha)*cos(alpha)
1908: ! py=cos(alpha)*sin(beta)1465: ! py=cos(alpha)*sin(beta)
1909: ! pz=sin(alpha)1466: ! pz=sin(alpha)
1910: !use commons1467: !use commons
1911: implicit none1468: implicit none
1912: 1469: 
1913: DOUBLE PRECISION        :: px1, py1, pz1,px,py,pz1470: DOUBLE PRECISION        :: px1, py1, pz1,px,py,pz
1914: 1471: 


r31914/key.f90 2017-02-15 10:30:13.980849470 +0000 r31913/key.f90 2017-02-15 10:30:14.872861366 +0000
 49:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, & 49:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, &
 50:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, & 50:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, &
 51:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 51:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 52:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 52:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 53:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 53:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 54:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 54:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 55:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 55:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 56:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 56:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 57:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, & 57:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, &
 58:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, & 58:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, &
 59:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T 59:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT
 60:  60: 
 61: ! sy349 > for testing the flatpath after dneb 61: ! sy349 > for testing the flatpath after dneb
 62:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:) 62:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:)
 63:       LOGICAL FLATPATHT 63:       LOGICAL FLATPATHT
 64:  64: 
 65: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 65: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 66:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 66:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 67:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 67:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 68:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 68:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
 69:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 69:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
108:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &108:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &
109:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &109:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &
110:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2110:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2
111: 111: 
112: !     sf344112: !     sf344
113:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &113:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &
114:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &114:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &
115:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT115:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT
116:  116:  
117:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)117:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)
118:       DOUBLE PRECISION, ALLOCATABLE :: PYADDREP(:,:), PYADDATT(:,:) 
119:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)118:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)
120:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:),SITE(:,:)119:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:),SITE(:,:)
121:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF120:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF
122:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET,PYT,PYCOLDFUSION121:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET,PYT,PYCOLDFUSION
123:       INTEGER          :: MAXINTERACTIONS,PYBINARYTYPE1,MYUNIT,NPYSITE122:       INTEGER          :: MAXINTERACTIONS,PYBINARYTYPE1,MYUNIT,NPYSITE
124:       DOUBLE PRECISION :: GBANISOTROPYR, GBWELLDEPTHR123:       DOUBLE PRECISION :: GBANISOTROPYR, GBWELLDEPTHR
125: !     CHRIS124: !     CHRIS
126: !     MCP125: !     MCP
127:       LOGICAL :: AMHT126:       LOGICAL :: AMHT
128:       CHARACTER(LEN=3) :: RES_TYPE_ALA='ALA'127:       CHARACTER(LEN=3) :: RES_TYPE_ALA='ALA'


r31914/keywords.f 2017-02-15 10:30:14.204852455 +0000 r31913/keywords.f 2017-02-15 10:30:15.096864354 +0000
730:          MIEF_CUTT=.FALSE.730:          MIEF_CUTT=.FALSE.
731:          MIEF_BOX(1:3) = 1.0D9731:          MIEF_BOX(1:3) = 1.0D9
732:          MIEF_RCUT= 1.0D9732:          MIEF_RCUT= 1.0D9
733:          !733:          !
734:          ! UNDOCUMENTED keywords/parameters734:          ! UNDOCUMENTED keywords/parameters
735:          ! 735:          ! 
736:          TWISTT=.FALSE.736:          TWISTT=.FALSE.
737:          LJADDT=.FALSE.737:          LJADDT=.FALSE.
738:          LJADD2T=.FALSE.738:          LJADD2T=.FALSE.
739:          PYADDT=.FALSE.739:          PYADDT=.FALSE.
740:          PYADD2T=.FALSE. 
741:          INVERTPT=.FALSE.740:          INVERTPT=.FALSE.
742:          DNEBEFRAC=0.0D0741:          DNEBEFRAC=0.0D0
743:          MORPHT=.FALSE.742:          MORPHT=.FALSE.
744:          GREATCIRCLET=.FALSE.743:          GREATCIRCLET=.FALSE.
745:          MAXTSENERGY=1.0D100744:          MAXTSENERGY=1.0D100
746:          MAXBARRIER=1.0D100745:          MAXBARRIER=1.0D100
747:          MAXMAXBARRIER=1.0D100746:          MAXMAXBARRIER=1.0D100
748:          ReoptimiseEndpoints=.False.747:          ReoptimiseEndpoints=.False.
749:          ANGLEAXIS=.FALSE.748:          ANGLEAXIS=.FALSE.
750:          NFAILMAX=2749:          NFAILMAX=2
5520:          CALL READF(PYSIGNOT)5519:          CALL READF(PYSIGNOT)
5521:          CALL READF(PYEPSNOT)5520:          CALL READF(PYEPSNOT)
5522: 5521: 
5523:          RBAAT = .TRUE.5522:          RBAAT = .TRUE.
5524:          ! Rigid body SITE, NRBSITES, NTSITES information specified in py_input routine5523:          ! Rigid body SITE, NRBSITES, NTSITES information specified in py_input routine
5525: 5524: 
5526:          ! Specify cutoff for potential in absolute units5525:          ! Specify cutoff for potential in absolute units
5527:          IF (NITEMS.GT.3) THEN5526:          IF (NITEMS.GT.3) THEN
5528:             CALL READF(PCUTOFF)5527:             CALL READF(PCUTOFF)
5529:             PARAMONOVCUTOFF=.TRUE.5528:             PARAMONOVCUTOFF=.TRUE.
5530:             WRITE(*,*) "multisitepy cutoff: ", PCUTOFF5529:             WRITE(MYUNIT,*) "multisitepy cutoff: ", PCUTOFF
5531:          ENDIF5530:          ENDIF
5532:          ! Specify periodic boundary conditions (PBCs)5531:          ! Specify periodic boundary conditions (PBCs)
5533:          IF (NITEMS.GT.4) THEN5532:          IF (NITEMS.GT.4) THEN
5534:             ! control which dimensions have periodic boundaries with a string 'XYZ', always put x before y before z.            5533:             ! control which dimensions have periodic boundaries with a string 'XYZ', always put x before y before z.            
5535:             ! eg ...  Xz 20 30  specifies PBC on X and Z directions.  The X box size will be 20, the Z box size 305534:             ! eg ...  Xz 20 30  specifies PBC on X and Z directions.  The X box size will be 20, the Z box size 30
5536:             CALL READA(PBC)5535:             CALL READA(PBC)
5537:             BOXLX=05536:             BOXLX=0
5538:             BOXLY=05537:             BOXLY=0
5539:             BOXLZ=05538:             BOXLZ=0
5540:             IF (SCAN(PBC,'Xx').NE.0) THEN5539:             IF (SCAN(PBC,'Xx').NE.0) THEN
5541:                 PARAMONOVPBCX=.TRUE.5540:                 PARAMONOVPBCX=.TRUE.
5542:                 CALL READF(BOXLX)5541:                 CALL READF(BOXLX)
5543:                 BOXLX = BOXLX*PCUTOFF5542:                 BOXLX = BOXLX*PCUTOFF
5544:                 WRITE(*,*) "PBC X:",BOXLX5543:                 WRITE(MYUNIT,*) "PBC X:",BOXLX
5545:             ENDIF5544:             ENDIF
5546:             IF (SCAN(PBC,'Yy').NE.0) THEN5545:             IF (SCAN(PBC,'Yy').NE.0) THEN
5547:                 PARAMONOVPBCY=.TRUE.5546:                 PARAMONOVPBCY=.TRUE.
5548:                 CALL READF(BOXLY)5547:                 CALL READF(BOXLY)
5549:                 BOXLY = BOXLY*PCUTOFF5548:                 BOXLY = BOXLY*PCUTOFF
5550:                 WRITE(*,*) "PBC Y:",BOXLY5549:                 WRITE(MYUNIT,*) "PBC Y:",BOXLY
5551:             ENDIF5550:             ENDIF
5552:             IF (SCAN(PBC,'Zz').NE.0) THEN5551:             IF (SCAN(PBC,'Zz').NE.0) THEN
5553:                 PARAMONOVPBCZ=.TRUE.5552:                 PARAMONOVPBCZ=.TRUE.
5554:                 CALL READF(BOXLZ)5553:                 CALL READF(BOXLZ)
5555:                 BOXLZ = BOXLZ*PCUTOFF5554:                 BOXLZ = BOXLZ*PCUTOFF
5556:                 WRITE(*,*) "PBC Z:",BOXLZ5555:                 WRITE(MYUNIT,*) "PBC Z:",BOXLZ
5557:             ENDIF5556:             ENDIF
5558:          ENDIF5557:          ENDIF
5559:       ELSE IF (WORD .EQ. 'PYGRAVITY') THEN5558:       ELSE IF (WORD .EQ. 'PYGRAVITY') THEN
5560:            CALL READF(PYGRAVITYC1)5559:            CALL READF(PYGRAVITYC1)
5561:            CALL READF(PYGRAVITYC2)5560:            CALL READF(PYGRAVITYC2)
5562:            EFIELDT=.TRUE.5561:            EFIELDT=.TRUE.
5563:       ELSE IF (WORD.EQ.'PYOVERLAPTHRESH') THEN5562:       ELSE IF (WORD.EQ.'PYOVERLAPTHRESH') THEN
5564:             CALL READF(PYOVERLAPTHRESH)5563:             CALL READF(PYOVERLAPTHRESH)
5565:             PYLOCALSTEP(:)=1.0D05564:             PYLOCALSTEP(:)=1.0D0
5566:             WRITE(*,'(A,F8.3)') 'keywords> ellipsoids considered to overlap for an ECF value of ', PYOVERLAPTHRESH5565:             WRITE(*,'(A,F8.3)') 'keywords> ellipsoids considered to overlap for an ECF value of ', PYOVERLAPTHRESH
5580:          LUNIT=GETUNIT()5579:          LUNIT=GETUNIT()
5581:          OPEN(LUNIT,FILE='epsilon',STATUS='OLD')5580:          OPEN(LUNIT,FILE='epsilon',STATUS='OLD')
5582:          IF (.NOT.ALLOCATED(PYADDEPS)) ALLOCATE(PYADDEPS(NADDTARGET,NADDTARGET))5581:          IF (.NOT.ALLOCATED(PYADDEPS)) ALLOCATE(PYADDEPS(NADDTARGET,NADDTARGET))
5583:          DO J1=1,NADDTARGET5582:          DO J1=1,NADDTARGET
5584:             DO J2=1,NADDTARGET5583:             DO J2=1,NADDTARGET
5585:                READ(LUNIT,*) PYADDEPS(J2,J1)5584:                READ(LUNIT,*) PYADDEPS(J2,J1)
5586:                WRITE(*,'(2I6,G20.10)') J1,J2,PYADDEPS(J2,J1)5585:                WRITE(*,'(2I6,G20.10)') J1,J2,PYADDEPS(J2,J1)
5587:             ENDDO5586:             ENDDO
5588:          ENDDO5587:          ENDDO
5589:          CLOSE(LUNIT)5588:          CLOSE(LUNIT)
5590:       ELSE IF (WORD.EQ.'PYADD2') THEN 
5591:          PYADD2T=.TRUE. 
5592:          CALL READI(NADDTARGET) 
5593:          WRITE(*,'(A,I6)') 'keyword> Target cluster size is ',NADDTARGET 
5594:          IF (MOD(NATOMS/2,NADDTARGET).NE.0) THEN 
5595:             WRITE(*,'(A,I6)') 'keyword> ERROR, target cluster size is not a factor of the number of PY particles ', 
5596:      &          NATOMS/2 
5597:             STOP 
5598:          ENDIF 
5599:          LUNIT=GETUNIT() 
5600:          OPEN(LUNIT,FILE='epsilon',STATUS='OLD') 
5601:          IF (.NOT.ALLOCATED(PYADDREP)) ALLOCATE(PYADDREP(NADDTARGET,NADDTARGET)) 
5602:          IF (.NOT.ALLOCATED(PYADDATT)) ALLOCATE(PYADDATT(NADDTARGET,NADDTARGET)) 
5603:          DO J1=1,NADDTARGET 
5604:             DO J2=1,NADDTARGET 
5605:                READ(LUNIT,*) PYADDREP(J2,J1), PYADDATT(J2,J1) 
5606:                WRITE(*,'(2I6,2G20.10)') J1,J2,PYADDREP(J2,J1),PYADDATT(J2,J1) 
5607:             ENDDO 
5608:          ENDDO 
5609:          CLOSE(LUNIT) 
5610:  
5611:       ELSE IF (WORD.EQ.'PYGPERIODIC') THEN5589:       ELSE IF (WORD.EQ.'PYGPERIODIC') THEN
5612:             NRBSITES = 25590:             NRBSITES = 2
5613:             ALLOCATE(RBSITE(NRBSITES,3))5591:             ALLOCATE(RBSITE(NRBSITES,3))
5614:             PYGPERIODICT = .TRUE.5592:             PYGPERIODICT = .TRUE.
5615: !            ANGLEAXIS2=.TRUE.5593: !            ANGLEAXIS2=.TRUE.
5616:             RBAAT=.TRUE.5594:             RBAAT=.TRUE.
5617:             CALL READF(PYA1(1))5595:             CALL READF(PYA1(1))
5618:             CALL READF(PYA1(2))5596:             CALL READF(PYA1(2))
5619:             CALL READF(PYA1(3))5597:             CALL READF(PYA1(3))
5620:             CALL READF(PYA2(1))5598:             CALL READF(PYA2(1))


r31914/potential.f 2017-02-15 10:30:14.428855443 +0000 r31913/potential.f 2017-02-15 10:30:15.316867289 +0000
1817:             IF (SSTEST) THEN1817:             IF (SSTEST) THEN
1818:                CALL PYGSECDER(COORDS,SSTEST)1818:                CALL PYGSECDER(COORDS,SSTEST)
1819:             END IF1819:             END IF
1820:             IF (PTEST) THEN1820:             IF (PTEST) THEN
1821:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '1821:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
1822:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1822:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
1823:             END IF1823:             END IF
1824:          ELSE IF (PYGPERIODICT.OR.PYBINARYT) THEN1824:          ELSE IF (PYGPERIODICT.OR.PYBINARYT) THEN
1825:             IF (PYADDT) THEN1825:             IF (PYADDT) THEN
1826:                CALL PYGPERIODICADD(COORDS,VNEW,ENERGY,GTEST,SSTEST)1826:                CALL PYGPERIODICADD(COORDS,VNEW,ENERGY,GTEST,SSTEST)
1827:             ELSEIF (PYADD2T) THEN 
1828:                CALL PYGPERIODICADD2(COORDS,VNEW,ENERGY,GTEST,SSTEST) 
1829:             ELSE1827:             ELSE
1830:                CALL PYGPERIODIC (COORDS,VNEW,ENERGY,GTEST,SSTEST)1828:                CALL PYGPERIODIC (COORDS,VNEW,ENERGY,GTEST,SSTEST)
1831:             ENDIF1829:             ENDIF
1832:             ! CALL PARAMONOVNUMFIRSTDER (COORDS,VNEW,ENERGY,GTEST,SSTEST)1830:             ! CALL PARAMONOVNUMFIRSTDER (COORDS,VNEW,ENERGY,GTEST,SSTEST)
1833:             ! STOP1831:             ! STOP
1834:             ! WRITE(*,*) 'VNEW='1832:             ! WRITE(*,*) 'VNEW='
1835:             ! DO J1=1,3*NATOMS1833:             ! DO J1=1,3*NATOMS
1836:             ! WRITE(*,*) VNEW(J1)1834:             ! WRITE(*,*) VNEW(J1)
1837:             ! END DO1835:             ! END DO
1838:             ! STOP1836:             ! STOP


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0