hdiff output

r31852/efol.f90 2017-01-30 12:30:09.460022749 +0000 r31851/efol.f90 2017-01-30 12:30:13.684080189 +0000
371:             DO J1=1,NATOMS+3371:             DO J1=1,NATOMS+3
372:                ZT(J1)=.FALSE.372:                ZT(J1)=.FALSE.
373:             ENDDO373:             ENDDO
374:             NZERO=NATOMS+3374:             NZERO=NATOMS+3
375:          ENDIF375:          ENDIF
376:       ELSE IF (FREEZE) THEN376:       ELSE IF (FREEZE) THEN
377:          NZERO=3*NFREEZE377:          NZERO=3*NFREEZE
378:          DO J1=1,NZERO378:          DO J1=1,NZERO
379:             ZT(J1)=.FALSE.379:             ZT(J1)=.FALSE.
380:          ENDDO380:          ENDDO
381:       ELSEIF(GBT.OR.GBDT) THEN381:       ELSEIF (GBT .OR. GBDT .OR. (PYGT .AND. UNIAXT).OR.((PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT).OR.STOCKAAT) THEN
382: !.OR.(PYGT.AND.UNIAXT).OR.((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT).OR.STOCKAAT) THEN 
383: 382: 
384:          IF (EFIELDT) THEN383:          IF (EFIELDT) THEN
385:             NZERO = NATOMS/2 + 4384:             NZERO = NATOMS/2 + 4
386:          ELSE385:          ELSE
387:             NZERO = NATOMS/2 + 6 386:             NZERO = NATOMS/2 + 6 
388:          ENDIF387:          ENDIF
389:          DO J1=1,NZERO388:          DO J1=1,NZERO
390:             ZT(J1)=.FALSE.389:             ZT(J1)=.FALSE.
391:          ENDDO390:          ENDDO
392:       391:       


r31852/fetchz.f 2017-01-30 12:30:09.684025795 +0000 r31851/fetchz.f 2017-01-30 12:30:13.916083344 +0000
 84:  84: 
 85: C 85: C
 86: C  Thomson problem: 86: C  Thomson problem:
 87: C 87: C
 88:       DOUBLE PRECISION, PARAMETER :: HALFPI=1.570796327D0 88:       DOUBLE PRECISION, PARAMETER :: HALFPI=1.570796327D0
 89:       DOUBLE PRECISION THTEMP(3*NATOMS), DIST 89:       DOUBLE PRECISION THTEMP(3*NATOMS), DIST
 90:       DOUBLE PRECISION, ALLOCATABLE :: FOCK(:,:), MOCOEFF(:,:) 90:       DOUBLE PRECISION, ALLOCATABLE :: FOCK(:,:), MOCOEFF(:,:)
 91:  91: 
 92:       IF (CHRMMT.OR.UNRST.OR.AMBERT.OR.NABT.OR.AMBER12T.OR.OPEPT) THEN 92:       IF (CHRMMT.OR.UNRST.OR.AMBERT.OR.NABT.OR.AMBER12T.OR.OPEPT) THEN
 93:          IF (.NOT.ALLOCATED(IATNUM)) ALLOCATE (IATNUM(NATOMS))   ! ATMASS already set up 93:          IF (.NOT.ALLOCATED(IATNUM)) ALLOCATE (IATNUM(NATOMS))   ! ATMASS already set up
 94:       ELSE IF (PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.PYGT.OR.MULTISITEPYT) THEN 94:       ELSE IF (PYGPERIODICT.OR.PYBINARYT.OR.PYGT.OR.MULTISITEPYT) THEN
 95:          IF (.NOT.ALLOCATED(IATNUM)) ALLOCATE (IATNUM(NATOMS),ATMASS(NATOMS)) 95:          IF (.NOT.ALLOCATED(IATNUM)) ALLOCATE (IATNUM(NATOMS),ATMASS(NATOMS))
 96:          ATMASS(:)=1.0D0 96:          ATMASS(:)=1.0D0
 97:       ELSE 97:       ELSE
 98:          IF (.NOT.ALLOCATED(IATNUM)) ALLOCATE (IATNUM(NATOMS)) 98:          IF (.NOT.ALLOCATED(IATNUM)) ALLOCATE (IATNUM(NATOMS))
 99:          IF (.NOT.ALLOCATED(ATMASS)) ALLOCATE (ATMASS(NATOMS)) 99:          IF (.NOT.ALLOCATED(ATMASS)) ALLOCATE (ATMASS(NATOMS))
100: C        sf344> initialise ATMASS values to not default to zero 100: C        sf344> initialise ATMASS values to not default to zero 
101: C               uninitialised stuff (may cause problems in INERTIA otherwise)101: C               uninitialised stuff (may cause problems in INERTIA otherwise)
102:          ! With RIGIDINIT, masses should now be specified by a separate file and will be initialised102:          ! With RIGIDINIT, masses should now be specified by a separate file and will be initialised
103:          ! in GENRIGID_INITIALISE.103:          ! in GENRIGID_INITIALISE.
104:          IF (.NOT. RIGIDINIT) ATMASS(:)=1.0D0104:          IF (.NOT. RIGIDINIT) ATMASS(:)=1.0D0


r31852/gay-berne.f90 2017-01-30 12:30:09.912028895 +0000 r31851/gay-berne.f90 2017-01-30 12:30:14.148086498 +0000
141:            END IF141:            END IF
142:           END DO142:           END DO
143:         END DO143:         END DO
144: !        WRITE(*,*) 'epsilon1: ', epsilon1(:,:,:)144: !        WRITE(*,*) 'epsilon1: ', epsilon1(:,:,:)
145:        END IF145:        END IF
146: ! cutoff stuff for the extra LJ sites146: ! cutoff stuff for the extra LJ sites
147:        cut = PCUTOFF*PCUTOFF ! cutoff squared147:        cut = PCUTOFF*PCUTOFF ! cutoff squared
148:        ron = PCUTOFF*0.9D0148:        ron = PCUTOFF*0.9D0
149:        ron2 = ron*ron149:        ron2 = ron*ron
150:        range2inv3=1.0D0/(cut-ron2)**3150:        range2inv3=1.0D0/(cut-ron2)**3
151: !     FROM IUT PARAMETERS151: !     FROM INPUT PARAMETERS
152: 152: 
153:       REALNATOMS = NATOMS/2153:       REALNATOMS = NATOMS/2
154:       OFFSET     = 3*REALNATOMS154:       OFFSET     = 3*REALNATOMS
155:       155:       
156:       DO K1 = 1, 3156:       DO K1 = 1, 3
157:         I3(K1,K1) = 1.0D0157:         I3(K1,K1) = 1.0D0
158:       ENDDO158:       ENDDO
159:       DO J1=1,REALNATOMS159:       DO J1=1,REALNATOMS
160:        DO K1 = 1, 3160:        DO K1 = 1, 3
161:          AEZR1(J1,K1,K1) = 1.D0/(PYA1bin(J1,K1)*PYA1bin(J1,K1))161:          AEZR1(J1,K1,K1) = 1.D0/(PYA1bin(J1,K1)*PYA1bin(J1,K1))
1188: &                        ' with the number of atoms determined from the path.xyz file!!!'1188: &                        ' with the number of atoms determined from the path.xyz file!!!'
1189:         STOP1189:         STOP
1190:    END IF1190:    END IF
1191:    DO J2=1,REALNATOMS1191:    DO J2=1,REALNATOMS
1192:         READ(133,'(A1,6X,6G20.10)') ATOMLABELS(J2),COORDSNEW(3*J2-2),COORDSNEW(3*J2-1),COORDSNEW(3*J2), &1192:         READ(133,'(A1,6X,6G20.10)') ATOMLABELS(J2),COORDSNEW(3*J2-2),COORDSNEW(3*J2-1),COORDSNEW(3*J2), &
1193:         & COORDSNEW(3*(J2+REALNATOMS)-2),COORDSNEW(3*(J2+REALNATOMS)-1),COORDSNEW(3*(J2+REALNATOMS))1193:         & COORDSNEW(3*(J2+REALNATOMS)-2),COORDSNEW(3*(J2+REALNATOMS)-1),COORDSNEW(3*(J2+REALNATOMS))
1194:    END DO1194:    END DO
1195:    COORDSDUMMY(:)=COORDSNEW(:)1195:    COORDSDUMMY(:)=COORDSNEW(:)
1196: !   COORDSDUMMY(50)=COORDSNEW(50)*10.0D01196: !   COORDSDUMMY(50)=COORDSNEW(50)*10.0D0
1197:    CALL MINPERMDIST(COORDSREF,COORDSDUMMY,NATOMS,.false.,1.0,1.0,1.0,.false.,.false.,D,DIST2,.false.,RMAT)1197:    CALL MINPERMDIST(COORDSREF,COORDSDUMMY,NATOMS,.false.,1.0,1.0,1.0,.false.,.false.,D,DIST2,.false.,RMAT)
1198:    D=D**2 ! since MINPERMDIST now returns the distance1198:    D=D**2 ! since minpermdist now returns the distance
1199:    COORDSNEW(:)=COORDSDUMMY(:)1199:    COORDSNEW(:)=COORDSDUMMY(:)
1200:    WRITE(*,*) 'writing frame ', J1+11200:    WRITE(*,*) 'writing frame ', J1+1
1201:    WRITE(134,*) REALNATOMS1201:    WRITE(134,*) REALNATOMS
1202:    WRITE(134,*) dummyline1202:    WRITE(134,*) dummyline
1203:    DO J2=1,REALNATOMS1203:    DO J2=1,REALNATOMS
1204:      WRITE(134,'(A1,6X,6G20.10)') ATOMLABELS(J2),COORDSNEW(3*J2-2),COORDSNEW(3*J2-1),COORDSNEW(3*J2), &1204:      WRITE(134,'(A1,6X,6G20.10)') ATOMLABELS(J2),COORDSNEW(3*J2-2),COORDSNEW(3*J2-1),COORDSNEW(3*J2), &
1205:         & COORDSNEW(3*(J2+REALNATOMS)-2),COORDSNEW(3*(J2+REALNATOMS)-1),COORDSNEW(3*(J2+REALNATOMS))1205:         & COORDSNEW(3*(J2+REALNATOMS)-2),COORDSNEW(3*(J2+REALNATOMS)-1),COORDSNEW(3*(J2+REALNATOMS))
1206:    END DO1206:    END DO
1207:    COORDSREF(:)=COORDSNEW(:)1207:    COORDSREF(:)=COORDSNEW(:)
1208:  END DO1208:  END DO
1322:   READ(nunit,*,iostat=iostatus) check1322:   READ(nunit,*,iostat=iostatus) check
1323: !  write(*,*) check,nunit1323: !  write(*,*) check,nunit
1324: end do1324: end do
1325:   nlines = nlines - 11325:   nlines = nlines - 1
1326:   REWIND(nunit)1326:   REWIND(nunit)
1327: RETURN1327: RETURN
1328: 1328: 
1329: 1329: 
1330: END SUBROUTINE DETERMINELINES1330: END SUBROUTINE DETERMINELINES
1331: 1331: 
1332:       SUBROUTINE TAKESTEPELPSD (COORDS) 
1333:  
1334: !     THIS ROUTINE TAKES STEP FOR SINGLE-SITE ELLIPSOIDAL BODIES ENSURING NO OVERLAP 
1335:       USE commons, only : natoms 
1336:       USE key 
1337:       use pymodule, only : angle,angle2,twopi 
1338:  
1339:       IMPLICIT NONE 
1340:  
1341:       INTEGER          :: NP, JMAX, JMAX2, REALNATOMS, OFFSET 
1342:       INTEGER          :: J1, J2, J3, J4, J5, J6 
1343: !      INTEGER          :: PTINDX, I, J, K 
1344:       LOGICAL          :: OVRLPT 
1345:  
1346:       DOUBLE PRECISION :: PI, DUMMY, DUMMY2, SAVECOORDS(3*NATOMS), COORDS(3*NATOMS), COORDSO(3*NATOMS) 
1347:       DOUBLE PRECISION :: DIST(3*NATOMS/2), XMASS, YMASS, ZMASS, DMAX, VMAX, VMAX2, VAT(3*NATOMS),VATO(3*NATOMS) 
1348:       DOUBLE PRECISION :: VMIN, CMMAX, CMDIST(NATOMS/2), LOCALSTEP, STEP, ASTEP, OSTEP 
1349:       DOUBLE PRECISION :: DPRAND, RANDOM, THETA, PHI, BOXLX, BOXLY, BOXLZ 
1350:       DOUBLE PRECISION :: AE(3,3), BE(3,3)  
1351:       DOUBLE PRECISION :: RI(3), RJ(3), RIJ(3) 
1352:       DOUBLE PRECISION :: P1(3), P2(3), ABSRIJ, RCUT, XMIN, FMIN, ECFVAL 
1353: ! sf344 additions 
1354:       LOGICAL          :: DISSOCIATED(NATOMS/2),DISSOC,TMOVE,OMOVE,LJCAPSIDT 
1355:       INTEGER          :: CLOSESTATOMINDEX, K1, K2 
1356:       DOUBLE PRECISION :: MINDISTANCE, x1, x2, y1, y2, z1, z2, RCUTARRAY(NATOMS/2), RVEC(NATOMS/2,NATOMS/2) 
1357: !      DOUBLE PRECISION, ALLOCATABLE :: RVEC(:) 
1358:  
1359: !      WRITE(*,*) 'NATOMS in takestepelpsd ', NATOMS 
1360: !      IF(ALLOCATED(RVEC)) DEALLOCATE(RVEC) 
1361:  
1362:       PI         = 4.D0*ATAN(1.0D0) 
1363:       IF (GBT .OR. GBDT) THEN 
1364:  
1365:          RCUT       = GBEPSNOT*MAX(GBKAPPA,1.D0) 
1366:  
1367:       ELSE IF (PYGPERIODICT.OR.PYBINARYT) THEN 
1368:         DO J1=1,NATOMS/2 
1369:          RCUTARRAY(J1) = 2.0D0 * MAXVAL(PYA1bin(J1,:))  
1370:         END DO 
1371:          RCUT = MAXVAL(RCUTARRAY) 
1372:       ELSE IF (LJCAPSIDT) THEN 
1373:         PYA1bin(:,:)=0.5D0 
1374:         DO J1=1,NATOMS/2 
1375:          RCUTARRAY(J1) = 2.0D0 
1376:         END DO 
1377:          RCUT = 2.0D0 
1378:       ENDIF 
1379:  
1380:       REALNATOMS = NATOMS/2 
1381:       OFFSET     = 3 * REALNATOMS 
1382:  
1383:                   SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS) 
1384:  
1385:       DO J1 = 1,REALNATOMS 
1386:  
1387:          J2     = 3*J1 
1388:         IF (PARAMONOVPBCX) THEN 
1389: !ensure x component of particle 1 vector is within BoxLx/2 of zero. 
1390: !if it isn't then subtract integer number of boxlx's such that it is. 
1391:                 COORDS(J2-2)=COORDS(J2-2)-BOXLX*NINT(COORDS(J2-2)/BOXLX) 
1392:         ENDIF 
1393:  
1394:         IF (PARAMONOVPBCY) THEN 
1395: !ensure y component of particle 1 vector is within BoxLy/2 of zero. 
1396: !if it isn't then subtract integer number of boxly's such that it is. 
1397:                 COORDS(J2-1)=COORDS(J2-1)-BOXLY*NINT(COORDS(J2-1)/BOXLY) 
1398:         END IF 
1399:  
1400:         IF (PARAMONOVPBCZ) THEN 
1401: !ensure z component of particle 1 vector is within BoxLz/2 of zero. 
1402: !if it isn't then subtract integer number of boxlz's such that it is. 
1403:                 COORDS(J2)=COORDS(J2)-BOXLZ*NINT(COORDS(J2)/BOXLZ) 
1404:         ENDIF 
1405:  
1406:          DUMMY2 = COORDS(J2-2)**2 + COORDS(J2-1)**2 + COORDS(J2)**2 
1407:  
1408:          IF (DUMMY2 .GT. RADIUS) THEN 
1409: ! bring back the molecule within the radius 
1410:                 COORDS(J2-2)=COORDS(J2-2)-SQRT(RADIUS)*NINT(COORDS(J2-2)/SQRT(RADIUS)) 
1411:                 COORDS(J2-1)=COORDS(J2-1)-SQRT(RADIUS)*NINT(COORDS(J2-1)/SQRT(RADIUS)) 
1412:                 COORDS(J2)=COORDS(J2)-SQRT(RADIUS)*NINT(COORDS(J2)/SQRT(RADIUS)) 
1413:     WRITE(MYUNIT,'(A,2F20.10)') 'initial coordinate outside container -- bringing molecule back within the container radius' 
1414:  
1415: !            WRITE(*,'(A,I5,5F20.10)') 'J1,RAD,R**2,x,y,z:', J1, RADIUS, DUMMY2, COORDS(J2-2), & 
1416: !                                       COORDS(J2-1), COORDS(J2) 
1417: !            PRINT*, 'initial coordinate outside container -- increase container radius' 
1418: !            STOP 
1419:          END IF 
1420:  
1421:       END DO 
1422:  
1423:       DO J1 = 1,3*NATOMS 
1424:          COORDSO(J1) = COORDS(J1) 
1425:       END DO 
1426:  
1427:       DO J1 = 1,NATOMS/2 
1428:          VATO(J1) = VAT(J1) 
1429:       END DO 
1430: !     FIND THE CENTRE OF MASS 
1431:  
1432:       XMASS = 0.0D0 
1433:       YMASS = 0.0D0 
1434:       ZMASS = 0.0D0 
1435:  
1436:       DO J1 = 1,NATOMS/2 
1437:  
1438:          XMASS = XMASS + COORDS(3*(J1-1)+1) 
1439:          YMASS = YMASS + COORDS(3*(J1-1)+2) 
1440:          ZMASS = ZMASS + COORDS(3*(J1-1)+3) 
1441:  
1442:       ENDDO 
1443:  
1444:       XMASS = XMASS/(REALNATOMS) 
1445:       YMASS = YMASS/(REALNATOMS) 
1446:       ZMASS = ZMASS/(REALNATOMS) 
1447:  
1448: !     Find the most weakly bound atom, JMAX, the second most weakly bound atom, JMAX2, 
1449: !     and the pair energy of the most tightly bound atom, VMIN. An angular step is 
1450: !     taken for JMAX if its pair energy is > ASTEP*VMIN putting the atom at a radius of 
1451: !     DMAX (or CMMAX from CM of the cluster). 
1452:       JMAX  = 1 
1453:       DMAX  =  1.0D0 
1454:       VMAX  = -1.0D3 
1455:       VMAX2 = -1.0D3 
1456:       VMIN  =  0.0D0 
1457:       CMMAX =  1.0D0 
1458:  
1459:       DO J1 = 1, REALNATOMS 
1460:  
1461:          J2 = 3*J1 
1462:          DIST(J1)   = DSQRT( COORDS(J2-2)**2 + COORDS(J2-1)**2 + COORDS(J2)**2) 
1463:          CMDIST(J1) = DSQRT((COORDS(J2-2)-XMASS)**2+(COORDS(J2-1)-YMASS)**2+(COORDS(J2)-ZMASS)**2) 
1464:          IF (CMDIST(J1) .GT. CMMAX) CMMAX = CMDIST(J1) 
1465:          IF (DIST(J1) .GT. DMAX) DMAX = DIST(J1) 
1466:          IF (VAT(J1) .GT. VMAX) THEN 
1467:             VMAX = VAT(J1) 
1468:             JMAX = J1 
1469:          ELSE IF ((VAT(J1).LT. VMAX) .AND. (VAT(J1) .GT. VMAX2)) THEN 
1470:               VMAX2 = VAT(J1) 
1471:               JMAX2 = J1 
1472:          ENDIF 
1473:          IF (VAT(J1) .LT. VMIN) VMIN = VAT(J1) 
1474:  
1475:       ENDDO 
1476:  
1477:       IF (VAT(JMAX) > (ASTEP*VMIN) .AND. (.NOT.NORESET)) THEN 
1478:  
1479:          J2 = 3*JMAX 
1480:          THETA           = DPRAND()*PI 
1481:          PHI             = DPRAND()*PI*2.0D0 
1482:          COORDS(J2-2) = XMASS + (CMMAX+1.0D0)*DSIN(THETA)*DCOS(PHI) 
1483:          COORDS(J2-1) = YMASS + (CMMAX+1.0D0)*DSIN(THETA)*DSIN(PHI) 
1484:          COORDS(J2)   = ZMASS + (CMMAX+1.0D0)*DCOS(THETA) 
1485:          DUMMY           = COORDS(J2-2)**2 + COORDS(J2-1)**2 + COORDS(J2)**2 
1486:  
1487:          IF (DUMMY > RADIUS) THEN 
1488: !          RADIUS=1.0 
1489:             DUMMY           = DSQRT(RADIUS*0.99D0/DUMMY) 
1490:             COORDS(J2-2) = COORDS(J2-2)*DUMMY 
1491:             COORDS(J2-1) = COORDS(J2-1)*DUMMY 
1492:             COORDS(J2)   = COORDS(J2)*DUMMY 
1493:          END IF 
1494:  
1495:       ENDIF 
1496:  
1497:       DO J1 = 1, REALNATOMS 
1498:  
1499:          J3 = 3*J1 
1500:          J5 = OFFSET + J3 
1501:  
1502: !     CHECK FOR OVERLAP 
1503:  
1504:          OVRLPT = .TRUE. 
1505: !        WRITE(*,*) 'in here' 
1506:  
1507:  
1508:  
1509: 95       DO WHILE (OVRLPT) 
1510:             LOCALSTEP = 0.0D0 
1511:             IF (TMOVE) LOCALSTEP = STEP 
1512:  
1513: !            IF(FROZEN(J1))  LOCALSTEP = 0.0D0 
1514:  
1515:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1516:             COORDS(J3-2) = COORDS(J3-2) + LOCALSTEP*RANDOM 
1517:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1518:             COORDS(J3-1) = COORDS(J3-1) + LOCALSTEP*RANDOM 
1519:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1520:             COORDS(J3)   = COORDS(J3) + LOCALSTEP*RANDOM 
1521:  
1522:             LOCALSTEP = 0.0D0 
1523:             IF (OMOVE) LOCALSTEP = OSTEP  
1524:  
1525:             IF(PYBINARYT) LOCALSTEP = 0.0D0 
1526:  
1527:             IF(FROZEN(J1)) LOCALSTEP = 0.0D0 
1528:  
1529:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1530:             COORDS(J5-2) = COORDS(J5-2) + LOCALSTEP*RANDOM 
1531:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1532:             COORDS(J5-1) = COORDS(J5-1) + LOCALSTEP*RANDOM 
1533:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1534:             COORDS(J5)   = COORDS(J5) + LOCALSTEP*RANDOM 
1535: !    
1536:             OVRLPT = .FALSE. 
1537:  
1538:             RI(:)  = COORDS(J3-2:J3) 
1539:             P1(:)  = COORDS(J5-2:J5) 
1540:  
1541: !     ROTATION MATRIX 
1542:  
1543:             CALL ELPSDRTMT (P1, PYA1bin(J1,:), AE) 
1544:  
1545:             DO J2 = 1, REALNATOMS 
1546:  
1547:                IF (J2 == J1) CYCLE 
1548:   
1549:                J4    = 3*J2 
1550:                J6    = OFFSET + J4 
1551:                RJ(:) = COORDS(J4-2:J4) 
1552:                P2(:) = COORDS(J6-2:J6) 
1553:  
1554: !            IF(FROZEN(J2))  LOCALSTEP = 0.0D0 
1555:  
1556: !            RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1557: !            COORDS(J4-2) = COORDS(J4-2) + LOCALSTEP*RANDOM 
1558: !            RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1559: !            COORDS(J4-1) = COORDS(J4-1) + LOCALSTEP*RANDOM 
1560: !            RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1561: !            COORDS(J4)   = COORDS(J4) + LOCALSTEP*RANDOM 
1562:  
1563: !            LOCALSTEP = 0.0D0 
1564: !            IF (OMOVE) LOCALSTEP = OSTEP  
1565: !            IF(PYBINARYT) LOCALSTEP = 0.0D0 
1566:  
1567: !            IF(FROZEN(J2)) LOCALSTEP = 0.0D0 
1568:  
1569: !            RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1570: !            COORDS(J6-2) = COORDS(J6-2) + LOCALSTEP*RANDOM 
1571: !            RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1572: !            COORDS(J6-1) = COORDS(J6-1) + LOCALSTEP*RANDOM 
1573: !            RANDOM          = (DPRAND() - 0.5D0)*2.0D0 
1574: !            COORDS(J6)   = COORDS(J6) + LOCALSTEP*RANDOM 
1575:  
1576: !     ROTATION MATRIX 
1577:  
1578:                CALL ELPSDRTMT (P2, PYA1bin(J2,:), BE) 
1579:  
1580:                RIJ    = RI - RJ 
1581:                ABSRIJ = DSQRT(DOT_PRODUCT(RIJ,RIJ)) 
1582:                          
1583:                IF (ABSRIJ < RCUT) THEN  
1584:  
1585: !     DETERMINE ELLIPTIC CONTACT FUNCTION 
1586:  
1587:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE, BE, RIJ, XMIN, FMIN) 
1588:  
1589:                   ECFVAL = - FMIN 
1590: ! allow for a slight overlap  
1591:                   IF (ECFVAL < PYOVERLAPTHRESH) THEN 
1592: !                        WRITE(*,*) 'atoms overlapping', J1, J2, ECFVAL, ABSRIJ 
1593:                      OVRLPT = .TRUE. 
1594:                      GO TO 95 
1595:                   ENDIF 
1596:  
1597:                ENDIF 
1598:  
1599:             ENDDO  ! END LOOP OVER J2 
1600:          ENDDO  ! END WHILE 
1601:       ENDDO  ! END LOOP OVER J1    
1602:  
1603:                      IF (FREEZE) THEN 
1604:                         DO J2=1,NATOMS 
1605:                            IF (FROZEN(J2)) THEN 
1606:                               COORDS(3*(J2-1)+1:3*(J2-1)+3)=SAVECOORDS(3*(J2-1)+1:3*(J2-1)+3) 
1607:                            ENDIF 
1608:                         ENDDO 
1609:                      ENDIF 
1610:  
1611:  
1612: ! now determine whether any particle has dissociated 
1613:  IF(PYGPERIODICT) THEN 
1614:   DO 
1615:   DISSOC=.FALSE. 
1616:    DO J1=1,natoms/2-1 
1617:     J2=3*J1 
1618:     x1=COORDS(J2-2) 
1619:     y1=COORDS(J2-1) 
1620:     z1=COORDS(J2) 
1621:     
1622:     DISSOCIATED(J1)=.TRUE. 
1623:     DO K1=1,natoms/2 
1624:       K2=3*K1 
1625:       x2=COORDS(K2-2) 
1626:       y2=COORDS(K2-1) 
1627:       z2=COORDS(K2) 
1628:       RVEC(J1,K1)=SQRT((x2-x1)**2+(y2-y1)**2+(z2-z1)**2) 
1629: !      WRITE(*,*) J1,K1,RVEC(J1,K1), 4*RCUT 
1630:       IF(RVEC(J1,K1)<4*RCUT) DISSOCIATED(J1)=.FALSE. 
1631: !      WRITE(*,*) 'in here, cutoff=',pcutoff,J1 
1632:     END DO 
1633:     IF(DISSOCIATED(J1)) THEN 
1634:         ! determine the closest neighbor 
1635:         MINDISTANCE = HUGE(1.0D0) 
1636:         DO K1=1,natoms/2 
1637:          IF(K1==J1) CYCLE 
1638:          IF(RVEC(J1,K1)<MINDISTANCE) THEN 
1639:            MINDISTANCE=RVEC(J1,K1) 
1640:            CLOSESTATOMINDEX=K1 
1641:          END IF 
1642:         END DO 
1643:         WRITE(MYUNIT,'(A,I6,A,F10.6,I6)') 'Atom ', J1, ' is dissociated, distance to the closest atom: ', & 
1644:         & MINDISTANCE, CLOSESTATOMINDEX 
1645:         DISSOC=.TRUE. 
1646:         ! move the dissociated atom closer to the closest atom 
1647:         COORDS(J2-2)=0.5D0*(COORDS(3*CLOSESTATOMINDEX-2)+COORDS(J2-2)) 
1648:         COORDS(J2-1)=0.5D0*(COORDS(3*CLOSESTATOMINDEX-1)+COORDS(J2-1)) 
1649:         COORDS(J2  )=0.5D0*(COORDS(3*CLOSESTATOMINDEX  )+COORDS(J2  )) 
1650:     END IF 
1651:    END DO 
1652:    IF(.NOT.DISSOC) EXIT 
1653:   END DO 
1654:         DO J1=NATOMS/2+1,NATOMS 
1655:             J2=3*J1 
1656:             angle2=dot_product(COORDS(J2-2:J2),COORDS(J2-2:J2)) 
1657:  
1658:             if(angle2>twopi**2) then 
1659:                 angle2=sqrt(angle2) 
1660:                 angle = angle2 - int(angle2/twopi)*twopi 
1661:                 COORDS(J2-2:J2)=COORDS(J2-2:J2)/angle2 * angle 
1662:             end if 
1663:         END DO 
1664:   END IF 
1665: END SUBROUTINE TAKESTEPELPSD 
1666:  
1667:       SUBROUTINE ELPSDRTMT (P, ESA, A) 
1668:  
1669:       IMPLICIT NONE 
1670:  
1671:       INTEGER          :: K 
1672:       DOUBLE PRECISION :: E(3,3), I3(3,3), ROT2(3,3), ROTMAT(3,3), P(3), THETA, THETA2  
1673:       DOUBLE PRECISION :: AEZR1(3,3), A(3,3), ESA(3) 
1674:  
1675:       I3(:,:)    = 0.D0 
1676:       AEZR1(:,:) = 0.D0 
1677:       I3(1,1)    = 1.D0 
1678:       I3(2,2)    = 1.D0 
1679:       I3(3,3)    = 1.D0 
1680:         
1681:       THETA   = DSQRT(DOT_PRODUCT(P,P)) 
1682:  
1683:       IF (THETA == 0.D0) THEN 
1684:  
1685:          ROTMAT = I3 
1686:  
1687:       ELSE 
1688:  
1689:          THETA2  = THETA * THETA 
1690:          E(:,:)  = 0.D0 
1691:          E(1,2)  = -P(3) 
1692:          E(1,3)  =  P(2) 
1693:          E(2,3)  = -P(1) 
1694:          E(2,1)  = -E(1,2) 
1695:          E(3,1)  = -E(1,3) 
1696:          E(3,2)  = -E(2,3) 
1697:          E       = E/THETA 
1698:  
1699:          ROTMAT = I3 + (1.D0-COS(THETA))*MATMUL(E,E) + E*SIN(THETA) 
1700:  
1701:       ENDIF 
1702:  
1703:       DO K = 1, 3 
1704:          AEZR1(K,K) = 1.D0/(ESA(K)*ESA(K)) 
1705:       ENDDO 
1706:  
1707:       A = MATMUL(ROTMAT,(MATMUL(AEZR1,(TRANSPOSE(ROTMAT))))) 
1708:  
1709:       END SUBROUTINE ELPSDRTMT  


r31852/geopt.f 2017-01-30 12:30:10.156032213 +0000 r31851/geopt.f 2017-01-30 12:30:14.372089545 +0000
405:             IF (BULKT.AND.TWOD) NEXMODES=NATOMS+2405:             IF (BULKT.AND.TWOD) NEXMODES=NATOMS+2
406:             IF (PULLT.OR.EFIELDT) NEXMODES=4406:             IF (PULLT.OR.EFIELDT) NEXMODES=4
407: ! hk286407: ! hk286
408:             IF (GTHOMSONT .AND. (GTHOMMET .EQ. 5) ) NEXMODES = 3408:             IF (GTHOMSONT .AND. (GTHOMMET .EQ. 5) ) NEXMODES = 3
409:             IF (GTHOMSONT .AND. (GTHOMMET < 5) ) NEXMODES = 1409:             IF (GTHOMSONT .AND. (GTHOMMET < 5) ) NEXMODES = 1
410:             IF (TWOD) NEXMODES=NEXMODES+NATOMS410:             IF (TWOD) NEXMODES=NEXMODES+NATOMS
411:             IF (FREEZE) THEN411:             IF (FREEZE) THEN
412:                NEXMODES=3*NFREEZE412:                NEXMODES=3*NFREEZE
413:             ENDIF413:             ENDIF
414:             IF (RBAAT) THEN414:             IF (RBAAT) THEN
415:              IF(UNIAXT) THEN 
416:                IF (EFIELDT) THEN 
417:                   NEXMODES = NATOMS/2 + 4 
418:                ELSE 
419:                   NEXMODES = NATOMS/2 + 6 
420:                ENDIF 
421:              ELSE 
422:                IF (EFIELDT) THEN415:                IF (EFIELDT) THEN
423:                   NEXMODES = 4416:                   NEXMODES = 4
424:                ELSE417:                ELSE
425:                   NEXMODES = 6418:                   NEXMODES = 6
426:                ENDIF419:                ENDIF
427:              END IF 
428:                IF (STOCKAAT) NEXMODES = NEXMODES + NATOMS/2420:                IF (STOCKAAT) NEXMODES = NEXMODES + NATOMS/2
429:             ENDIF421:             ENDIF
430:             IF (RINGPOLYMERT) THEN422:             IF (RINGPOLYMERT) THEN
431:                IF (RPSYSTEM(1:4).EQ.'AECK') THEN423:                IF (RPSYSTEM(1:4).EQ.'AECK') THEN
432:                   NEXMODES=0424:                   NEXMODES=0
433:                ELSE425:                ELSE
434:                   NEXMODES=6426:                   NEXMODES=6
435:                ENDIF427:                ENDIF
436:                IF (BFGSTST.OR.(INR.EQ.2)) NEXMODES=NEXMODES+1428:                IF (BFGSTST.OR.(INR.EQ.2)) NEXMODES=NEXMODES+1
437:                IF (BFGSMINT) NEXMODES=NEXMODES+2429:                IF (BFGSMINT) NEXMODES=NEXMODES+2
438:             ENDIF430:             ENDIF
439:             IF (STEALTHYT) THEN431:             IF (STEALTHYT) THEN
440:                IF (ENERGY.LT.1.0D-8) THEN432:                IF (ENERGY.LT.1.0D-8) THEN
441:                   NEXMODES=NATOMS*3 - STM*2433:                   NEXMODES=NATOMS*3 - STM*2
442:                ELSE434:                ELSE
443:                   NEXMODES=3435:                   NEXMODES=3
444:                ENDIF436:                ENDIF
445:             ENDIF437:             ENDIF
446:             IF (VARIABLES.OR.MIEFT .OR. NOTRANSROTT) NEXMODES=0438:             IF (VARIABLES.OR.MIEFT .OR. NOTRANSROTT) NEXMODES=0
447:             IF (BFGSTST) NEXMODES=NEXMODES+1 !sf344> so that the log product for transition state geometries can be calculated and dumped correctly 
448:             WRITE(*,'(A,I6)') ' geopt> Number of zero/imaginary eigenvalues assumed to be ',NEXMODES439:             WRITE(*,'(A,I6)') ' geopt> Number of zero/imaginary eigenvalues assumed to be ',NEXMODES
449: 440: 
450: ! Keyword-specific block441: ! Keyword-specific block
451:             IF (LOWESTFRQT) THEN442:             IF (LOWESTFRQT) THEN
452: C443: C
453: C  Calculate lowest non-zero eigenvalue and dump to min.data.info file444: C  Calculate lowest non-zero eigenvalue and dump to min.data.info file
454: C  No mass-weighting here!445: C  No mass-weighting here!
455: C  'U' specifies that the upper triangle contains the Hessian.446: C  'U' specifies that the upper triangle contains the Hessian.
456: C  Need to save HESS before this call and restore afterwards.447: C  Need to save HESS before this call and restore afterwards.
457: C448: C


r31852/key.f90 2017-01-30 12:30:10.392035423 +0000 r31851/key.f90 2017-01-30 12:30:14.596092592 +0000
105:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, &105:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, &
106:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &106:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &
107:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &107:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &
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)
116:  116:  
117:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)117:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)
118:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:),SITE(:,:)118:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:)
119:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF119:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF
120:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET,PYT,PYCOLDFUSION120:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET
121:       INTEGER          :: MAXINTERACTIONS,PYBINARYTYPE1,MYUNIT,NPYSITE121:       INTEGER          :: MAXINTERACTIONS,PYBINARYTYPE1,MYUNIT,NPYSITE
122:       DOUBLE PRECISION :: GBANISOTROPYR, GBWELLDEPTHR122:       DOUBLE PRECISION :: GBANISOTROPYR, GBWELLDEPTHR
123: !     CHRIS123: !     CHRIS
124: !     MCP124: !     MCP
125:       LOGICAL :: AMHT125:       LOGICAL :: AMHT
126:       CHARACTER(LEN=3) :: RES_TYPE_ALA='ALA'126:       CHARACTER(LEN=3) :: RES_TYPE_ALA='ALA'
127:       CHARACTER(LEN=5) :: TARFL127:       CHARACTER(LEN=5) :: TARFL
128:       CHARACTER(LEN=120) :: CASTEPJOB, CP2KJOB, ONETEPJOB, QCHEMJOB, VASPJOB, MULTISTART, MULTIFINISH,MOLPROJOB,REAXFFJOB128:       CHARACTER(LEN=120) :: CASTEPJOB, CP2KJOB, ONETEPJOB, QCHEMJOB, VASPJOB, MULTISTART, MULTIFINISH,MOLPROJOB,REAXFFJOB
129:       CHARACTER(LEN=30) :: QCHEMJOBPARAMS=''129:       CHARACTER(LEN=30) :: QCHEMJOBPARAMS=''
130:       CHARACTER(LEN=30) :: MOLPROJOBPARAMS=''130:       CHARACTER(LEN=30) :: MOLPROJOBPARAMS=''


r31852/keywords.f 2017-01-30 12:30:10.628038632 +0000 r31851/keywords.f 2017-01-30 12:30:14.832095801 +0000
875:          EYTRAPT=.FALSE.875:          EYTRAPT=.FALSE.
876:          MKTRAPT=.FALSE.876:          MKTRAPT=.FALSE.
877:          BFGSTSTOL=0.0001D0877:          BFGSTSTOL=0.0001D0
878:          EVPC=1.0D0878:          EVPC=1.0D0
879:          OHCELLT=.FALSE.879:          OHCELLT=.FALSE.
880:          CONDATT=.FALSE.880:          CONDATT=.FALSE.
881:          PAIRCOLOURT=.FALSE.881:          PAIRCOLOURT=.FALSE.
882:          NENDDUP=0882:          NENDDUP=0
883:          REVERSEUPHILLT=.FALSE.883:          REVERSEUPHILLT=.FALSE.
884:          ! sf344884:          ! sf344
885:          EFIELDT=.FALSE. 
886:          MACROIONT=.FALSE.885:          MACROIONT=.FALSE.
887:          NORMALMODET=.FALSE.886:          NORMALMODET=.FALSE.
888:          PYGPERIODICT=.FALSE.887:          PYGPERIODICT=.FALSE.
889:          PYT=.FALSE. 
890:          PYBINARYT=.FALSE.888:          PYBINARYT=.FALSE.
891:          MULTISITEPYT=.FALSE.889:          MULTISITEPYT=.FALSE.
892:          LJGSITET=.FALSE.890:          LJGSITET=.FALSE.
893:          LJSITE=.FALSE.891:          LJSITE=.FALSE.
894:          BLJSITE=.FALSE.892:          BLJSITE=.FALSE.
895:          LJSITECOORDST=.FALSE.893:          LJSITECOORDST=.FALSE.
896:          LJSITEATTR=.FALSE.894:          LJSITEATTR=.FALSE.
897:          PCUTOFF=999.0D0895:          PCUTOFF=999.0D0
898:          PARAMONOVCUTOFF=.FALSE. 
899:          CLOSESTALIGNMENT=.FALSE.896:          CLOSESTALIGNMENT=.FALSE.
900:          DF1T=.FALSE.897:          DF1T=.FALSE.
901:          PULLT=.FALSE.898:          PULLT=.FALSE.
902:          CHEMSHIFT=.FALSE.899:          CHEMSHIFT=.FALSE.
903:          METRICTENSOR=.FALSE.900:          METRICTENSOR=.FALSE.
904: 901: 
905:          FRQCONV = 1.0D0902:          FRQCONV = 1.0D0
906:          DUMMY_FRQCONV = 0.0D0903:          DUMMY_FRQCONV = 0.0D0
907: 904: 
908:          QUIPARGSTRT=.FALSE.905:          QUIPARGSTRT=.FALSE.
2809:             END IF2806:             END IF
2810:             IF(NITEMS.GT.4) THEN           ! binary ellipsoidal clusters will be set up only for two apex sites, not one2807:             IF(NITEMS.GT.4) THEN           ! binary ellipsoidal clusters will be set up only for two apex sites, not one
2811:                BLJSITE=.TRUE.               ! we also won't use the sigma parameter from now on, epsilon is enough for repulsive sites2808:                BLJSITE=.TRUE.               ! we also won't use the sigma parameter from now on, epsilon is enough for repulsive sites
2812:                CALL READF(PEPSILON1(2))2809:                CALL READF(PEPSILON1(2))
2813:                CALL READF(PSCALEFAC1(2))2810:                CALL READF(PSCALEFAC1(2))
2814:                CALL READF(PSCALEFAC2(2))2811:                CALL READF(PSCALEFAC2(2))
2815:                CALL READF(PEPSILON1(3))     ! this is epsilon for the interaction between A and B type ellipsoids2812:                CALL READF(PEPSILON1(3))     ! this is epsilon for the interaction between A and B type ellipsoids
2816:                MAXINTERACTIONS=3 ! attractive secondary apex sites not incorporated for binary systems2813:                MAXINTERACTIONS=3 ! attractive secondary apex sites not incorporated for binary systems
2817:                WRITE(*,'(A,3F8.3)') ' keyword> binary system with primary and secondary apex sites, ' //2814:                WRITE(*,'(A,3F8.3)') ' keyword> binary system with primary and secondary apex sites, ' //
2818:      &         'epsilon and heights for 2nd type particle: ', PEPSILON1(2), PSCALEFAC1(2), PSCALEFAC2(2)2815:      &         'epsilon and heights for 2nd type particle: ', PEPSILON1(2), PSCALEFAC1(2), PSCALEFAC2(2)
 2816: 
2819:             END IF2817:             END IF
2820:          ELSE IF (WORD.EQ.'EXTRALJSITEATTR') THEN2818:          ELSE IF (WORD.EQ.'EXTRALJSITEATTR') THEN
2821:             LJSITE=.TRUE.2819:             LJSITE=.TRUE.
2822:             LJSITEATTR=.TRUE.2820:             LJSITEATTR=.TRUE.
2823:             CALL READF(PSIGMAATTR(1))2821:             CALL READF(PSIGMAATTR(1))
2824:             CALL READF(PEPSILONATTR(1))2822:             CALL READF(PEPSILONATTR(1))
2825:             CALL READF(PSIGMAATTR(2))2823:             CALL READF(PSIGMAATTR(2))
2826:             CALL READF(PEPSILONATTR(2))2824:             CALL READF(PEPSILONATTR(2))
 2825: 
2827:             WRITE(*,'(A,4F8.3)') 'keyword> primary and secondary apex sites '//2826:             WRITE(*,'(A,4F8.3)') 'keyword> primary and secondary apex sites '//
2828:      &      'with normal LJ attraction, sigmas and epsilons: ',2827:      &      'with normal LJ attraction, sigmas and epsilons: ',
2829:      &      PSIGMAATTR(1), PEPSILONATTR(1), PSIGMAATTR(2), PEPSILONATTR(2)2828:      &      PSIGMAATTR(1), PEPSILONATTR(1), PSIGMAATTR(2), PEPSILONATTR(2)
2830:             MAXINTERACTIONS=42829:             MAXINTERACTIONS=4
2831:          ELSE IF (WORD.EQ.'LJSITECOORDS') THEN2830:          ELSE IF (WORD.EQ.'LJSITECOORDS') THEN
2832:             LJSITECOORDST=.TRUE.2831:             LJSITECOORDST=.TRUE.
2833:             CALL READF(LJSITECOORDS(1))2832:             CALL READF(LJSITECOORDS(1))
2834:             CALL READF(LJSITECOORDS(2))2833:             CALL READF(LJSITECOORDS(2))
2835:             CALL READF(LJSITECOORDS(3))2834:             CALL READF(LJSITECOORDS(3))
2836:          ELSE IF (WORD.EQ.'PYBINARY') THEN2835:          ELSE IF (WORD.EQ.'PYBINARY') THEN
5469:             CALL READF(PYSIGNOT)5468:             CALL READF(PYSIGNOT)
5470:             CALL READF(PYEPSNOT)5469:             CALL READF(PYEPSNOT)
5471: 5470: 
5472:             IF (PYA1(1) == PYA2(1) .AND. PYA1(2) == PYA2(2) .AND. PYA1(3) == PYA2(3)) THEN5471:             IF (PYA1(1) == PYA2(1) .AND. PYA1(2) == PYA2(2) .AND. PYA1(3) == PYA2(3)) THEN
5473:                RADIFT = .FALSE.5472:                RADIFT = .FALSE.
5474:             ELSE5473:             ELSE
5475:                RADIFT = .TRUE.5474:                RADIFT = .TRUE.
5476:             ENDIF5475:             ENDIF
5477: ! sf344> PY potential and extra LJ site5476: ! sf344> PY potential and extra LJ site
5478:          ELSE IF (WORD.EQ.'PYBINARY') THEN5477:          ELSE IF (WORD.EQ.'PYBINARY') THEN
5479:             NRBSITES = 25478:             NRBSITES = 1
5480:             ALLOCATE(RBSITE(NRBSITES,3))5479:             ALLOCATE(RBSITE(NRBSITES,3))
5481:             PYBINARYT=.TRUE.5480:             PYBINARYT=.TRUE.
5482:             ANGLEAXIS2=.TRUE.5481:             ANGLEAXIS2=.TRUE.
5483:             RBAAT=.TRUE.5482:             RBAAT=.TRUE.
5484:             RADIFT=.TRUE.5483:             RADIFT=.TRUE.
5485:             CALL READI(PYBINARYTYPE1)5484:             CALL READI(PYBINARYTYPE1)
5486:             CALL READF(PYA11(1))5485:             CALL READF(PYA11(1))
5487:             CALL READF(PYA11(2))5486:             CALL READF(PYA11(2))
5488:             CALL READF(PYA11(3))5487:             CALL READF(PYA11(3))
5489:             CALL READF(PYA21(1))5488:             CALL READF(PYA21(1))
5490:             CALL READF(PYA21(2))5489:             CALL READF(PYA21(2))
5491:             CALL READF(PYA21(3))5490:             CALL READF(PYA21(3))
5492:             CALL READF(PYA12(1))5491:             CALL READF(PYA12(1))
5493:             CALL READF(PYA12(2))5492:             CALL READF(PYA12(2))
5494:             CALL READF(PYA12(3))5493:             CALL READF(PYA12(3))
5495:             CALL READF(PYA22(1))5494:             CALL READF(PYA22(1))
5496:             CALL READF(PYA22(2))5495:             CALL READF(PYA22(2))
5497:             CALL READF(PYA22(3))5496:             CALL READF(PYA22(3))
5498:             CALL READF(PYSIGNOT)5497:             CALL READF(PYSIGNOT)
5499:             CALL READF(PYEPSNOT)5498:             CALL READF(PYEPSNOT)
 5499: 
5500:             IF(.NOT.ALLOCATED(PYA1bin)) ALLOCATE(PYA1bin(NATOMS/2,3))5500:             IF(.NOT.ALLOCATED(PYA1bin)) ALLOCATE(PYA1bin(NATOMS/2,3))
5501:             IF(.NOT.ALLOCATED(PYA2bin)) ALLOCATE(PYA2bin(NATOMS/2,3))5501:             IF(.NOT.ALLOCATED(PYA2bin)) ALLOCATE(PYA2bin(NATOMS/2,3))
5502:             DO J1=1,NATOMS/25502:             DO J1=1,NATOMS/2
5503:                IF(J1<=PYBINARYTYPE1) THEN5503:                IF(J1<=PYBINARYTYPE1) THEN
5504:                   PYA1bin(J1,:)=PYA11(:)5504:                   PYA1bin(J1,:)=PYA11(:)
5505:                   PYA2bin(J1,:)=PYA21(:)5505:                   PYA2bin(J1,:)=PYA21(:)
5506:                ELSE5506:                ELSE
5507:                   PYA1bin(J1,:)=PYA12(:)5507:                   PYA1bin(J1,:)=PYA12(:)
5508:                   PYA2bin(J1,:)=PYA22(:)5508:                   PYA2bin(J1,:)=PYA22(:)
5509:                END IF5509:                END IF
5510:             END DO5510:             END DO
5511:             RBSITE(1,1)=PYA11(1)5511:          ELSE IF (WORD.EQ.'PYOVERLAPTHRESH') THEN
5512:             RBSITE(2,1)=-PYA11(1) 
5513:             RBSITE(:,2:3)=0.0D0 
5514:       ELSE IF (WORD .EQ. 'PY') THEN 
5515:          ! Syntax: PY sig_0 eps_0 [cut] [XYZ boxx boxy boxz] 
5516:  
5517:          PYT = .TRUE. 
5518:          CALL READF(PYSIGNOT) 
5519:          CALL READF(PYEPSNOT) 
5520:  
5521:          RBAAT = .TRUE. 
5522:          ! Rigid body SITE, NRBSITES, NTSITES information specified in py_input routine 
5523:  
5524:          ! Specify cutoff for potential in absolute units 
5525:          IF (NITEMS.GT.3) THEN 
5526:             CALL READF(PCUTOFF) 
5527:             PARAMONOVCUTOFF=.TRUE. 
5528:             WRITE(MYUNIT,*) "multisitepy cutoff: ", PCUTOFF 
5529:          ENDIF 
5530:          ! Specify periodic boundary conditions (PBCs) 
5531:          IF (NITEMS.GT.4) THEN 
5532:             ! control which dimensions have periodic boundaries with a string 'XYZ', always put x before y before z.             
5533:             ! eg ...  Xz 20 30  specifies PBC on X and Z directions.  The X box size will be 20, the Z box size 30 
5534:             CALL READA(PBC) 
5535:             BOXLX=0 
5536:             BOXLY=0 
5537:             BOXLZ=0 
5538:             IF (SCAN(PBC,'Xx').NE.0) THEN 
5539:                 PARAMONOVPBCX=.TRUE. 
5540:                 CALL READF(BOXLX) 
5541:                 BOXLX = BOXLX*PCUTOFF 
5542:                 WRITE(MYUNIT,*) "PBC X:",BOXLX 
5543:             ENDIF 
5544:             IF (SCAN(PBC,'Yy').NE.0) THEN 
5545:                 PARAMONOVPBCY=.TRUE. 
5546:                 CALL READF(BOXLY) 
5547:                 BOXLY = BOXLY*PCUTOFF 
5548:                 WRITE(MYUNIT,*) "PBC Y:",BOXLY 
5549:             ENDIF 
5550:             IF (SCAN(PBC,'Zz').NE.0) THEN 
5551:                 PARAMONOVPBCZ=.TRUE. 
5552:                 CALL READF(BOXLZ) 
5553:                 BOXLZ = BOXLZ*PCUTOFF 
5554:                 WRITE(MYUNIT,*) "PBC Z:",BOXLZ 
5555:             ENDIF 
5556:          ENDIF 
5557:       ELSE IF (WORD .EQ. 'PYGRAVITY') THEN 
5558:            CALL READF(PYGRAVITYC1) 
5559:            CALL READF(PYGRAVITYC2) 
5560:            EFIELDT=.TRUE. 
5561:       ELSE IF (WORD.EQ.'PYOVERLAPTHRESH') THEN 
5562:             CALL READF(PYOVERLAPTHRESH)5512:             CALL READF(PYOVERLAPTHRESH)
5563:             PYLOCALSTEP(:)=1.0D05513:             PYLOCALSTEP(:)=1.0D0
5564:             WRITE(*,'(A,F8.3)') 'keywords> ellipsoids considered to overlap for an ECF value of ', PYOVERLAPTHRESH5514:             WRITE(*,'(A,F8.3)') 'keywords> ellipsoids considered to overlap for an ECF value of ', PYOVERLAPTHRESH
5565:             IF(NITEMS.GT.2) THEN5515:             IF(NITEMS.GT.2) THEN
5566:                CALL READF(PYLOCALSTEP(1))5516:                CALL READF(PYLOCALSTEP(1))
5567:                CALL READF(PYLOCALSTEP(2))5517:                CALL READF(PYLOCALSTEP(2))
5568:             END IF5518:             END IF
5569:       ELSE IF (WORD.EQ.'PYGPERIODIC') THEN5519:          ELSE IF (WORD.EQ.'PYGPERIODIC') THEN
5570:             NRBSITES = 2 
5571:             ALLOCATE(RBSITE(NRBSITES,3)) 
5572:             PYGPERIODICT = .TRUE.5520:             PYGPERIODICT = .TRUE.
5573: !            ANGLEAXIS2=.TRUE.5521:             ANGLEAXIS2=.TRUE.
5574:             RBAAT=.TRUE.5522:             RBAAT=.TRUE.
5575:             CALL READF(PYA1(1))5523:             CALL READF(PYA1(1))
5576:             CALL READF(PYA1(2))5524:             CALL READF(PYA1(2))
5577:             CALL READF(PYA1(3))5525:             CALL READF(PYA1(3))
5578:             CALL READF(PYA2(1))5526:             CALL READF(PYA2(1))
5579:             CALL READF(PYA2(2))5527:             CALL READF(PYA2(2))
5580:             CALL READF(PYA2(3))5528:             CALL READF(PYA2(3))
5581:             CALL READF(PYSIGNOT)5529:             CALL READF(PYSIGNOT)
5582:             CALL READF(PYEPSNOT)5530:             CALL READF(PYEPSNOT)
5583:             RBSITE(1,1)=PYA1(1)5531: 
5584:             RBSITE(2,1)=0.0D0 
5585:             RBSITE(:,2:3)=0.0D0 
5586:             NTSITES=NATOMS*NRBSITES/2 
5587:             IF(.NOT.ALLOCATED(PYA1bin)) ALLOCATE(PYA1bin(NATOMS/2,3))5532:             IF(.NOT.ALLOCATED(PYA1bin)) ALLOCATE(PYA1bin(NATOMS/2,3))
5588:             IF(.NOT.ALLOCATED(PYA2bin)) ALLOCATE(PYA2bin(NATOMS/2,3))5533:             IF(.NOT.ALLOCATED(PYA2bin)) ALLOCATE(PYA2bin(NATOMS/2,3))
5589:             DO J1=1,NATOMS/25534:             DO J1=1,NATOMS/2
5590:                PYA1bin(J1,:)=PYA1(:)5535:                PYA1bin(J1,:)=PYA1(:)
5591:                PYA2bin(J1,:)=PYA2(:)5536:                PYA2bin(J1,:)=PYA2(:)
5592:             END DO5537:             END DO
5593:             IF (PYA1(1) == PYA2(1) .AND. PYA1(2) == PYA2(2) .AND. PYA1(3) == PYA2(3)) THEN5538:             IF (PYA1(1) == PYA2(1) .AND. PYA1(2) == PYA2(2) .AND. PYA1(3) == PYA2(3)) THEN
5594:                RADIFT = .FALSE.5539:                RADIFT = .FALSE.
5595:             ELSE5540:             ELSE
5596:                RADIFT = .TRUE.5541:                RADIFT = .TRUE.


r31852/minpermdistrbcom.f90 2017-01-30 12:30:10.856041733 +0000 r31851/minpermdistrbcom.f90 2017-01-30 12:30:15.084099227 +0000
  1: !  1: !
  2: !     Follows the same procedure as in minpermdist.f90 to align a system of rigid bodies based  2: !     Follows the same procedure as in minpermdist.f90 to align a system of rigid bodies based
  3: !     on the centres of mass only.  3: !     on the centres of mass only.
  4: !     ----------------------------------------------------------------------------------------------  4: !     ----------------------------------------------------------------------------------------------
  5:    5:  
  6:       SUBROUTINE MINPERMDISTRBCOM(COORDSB,COORDSA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT)  6:       SUBROUTINE MINPERMDISTRBCOM(COORDSB,COORDSA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT)
  7: !     DISTANCE returns the squared distance  7: !     DISTANCE returns the squared distance
  8:   8: 
  9:       USE COMMONS, ONLY : NATOMS  9:       USE COMMONS, ONLY : NATOMS
 10:       USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, GEOMDIFFTOL, EFIELDT, NRBGROUP, PAPT, PTSTSTT, NOINVERSION, & 10:       USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, GEOMDIFFTOL, EFIELDT, NRBGROUP, PAPT, PTSTSTT, NOINVERSION 
 11: &                        MULTISITEPYT  
 12:  11: 
 13:       IMPLICIT NONE 12:       IMPLICIT NONE
 14:  13: 
 15:       INTEGER, PARAMETER :: MAXIMUMTRIES = 100 14:       INTEGER, PARAMETER :: MAXIMUMTRIES = 100
 16:       INTEGER            :: NPERM, PATOMS, NTRIES, NSIZE, JMAX, LOCMAX(1), J1, J2, J3, INFO 15:       INTEGER            :: NPERM, PATOMS, NTRIES, NSIZE, JMAX, LOCMAX(1), J1, J2, J3, INFO
 17:       INTEGER            :: INVERT, NORBIT1, NORBIT2, PERM(NATOMS), NCHOOSE2, NDUMMY, LPERM(NATOMS), NCHOOSE1 16:       INTEGER            :: INVERT, NORBIT1, NORBIT2, PERM(NATOMS), NCHOOSE2, NDUMMY, LPERM(NATOMS), NCHOOSE1
 18:       INTEGER            :: NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS) 17:       INTEGER            :: NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS)
 19:       DOUBLE PRECISION   :: COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DISTWP, DIST2, TEMPA(9*NATOMS)  18:       DOUBLE PRECISION   :: COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DISTWP, DIST2, TEMPA(9*NATOMS) 
 20:       DOUBLE PRECISION   :: DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS), DUMMYWP(3*NATOMS) 19:       DOUBLE PRECISION   :: DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS), DUMMYWP(3*NATOMS)
 21: !      DOUBLE PRECISION   :: XA(3*NATOMS*NRBSITES/2),  XB(3*NATOMS*NRBSITES/2), XBS(3*NATOMS*NRBSITES/2) 20: !      DOUBLE PRECISION   :: XA(3*NATOMS*NRBSITES/2),  XB(3*NATOMS*NRBSITES/2), XBS(3*NATOMS*NRBSITES/2)
100:          CMX = CMX + DUMMYA(J2-2) 99:          CMX = CMX + DUMMYA(J2-2)
101:          CMY = CMY + DUMMYA(J2-1)100:          CMY = CMY + DUMMYA(J2-1)
102:          CMZ = CMZ + DUMMYA(J2)101:          CMZ = CMZ + DUMMYA(J2)
103:       ENDDO102:       ENDDO
104: 103: 
105:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS104:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS
106:       DO J1 = 1, NATOMS105:       DO J1 = 1, NATOMS
107:          J2 = 3*J1106:          J2 = 3*J1
108:          DUMMYA(J2-2) = DUMMYA(J2-2) - CMX107:          DUMMYA(J2-2) = DUMMYA(J2-2) - CMX
109:          DUMMYA(J2-1) = DUMMYA(J2-1) - CMY108:          DUMMYA(J2-1) = DUMMYA(J2-1) - CMY
110: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY) 
111: !         IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)  
112:          DUMMYA(J2)   = DUMMYA(J2)   - CMZ109:          DUMMYA(J2)   = DUMMYA(J2)   - CMZ
113:       ENDDO110:       ENDDO
114: 111: 
115:       IF (EFIELDT .AND. INVERT ==-1) THEN112:       IF (EFIELDT .AND. INVERT ==-1) THEN
116:          RMATI(:,:) = 0.D0113:          RMATI(:,:) = 0.D0
117:          RMATI(1,1) = 1.D0; RMATI(2,2) =-1.D0; RMATI(3,3) = 1.D0114:          RMATI(1,1) = 1.D0; RMATI(2,2) =-1.D0; RMATI(3,3) = 1.D0
118:          DO J1 = 1, NATOMS115:          DO J1 = 1, NATOMS
119:             J2 = 3*J1116:             J2 = 3*J1
120:             DUMMYC(J2-2:J2) = MATMUL(RMATI,DUMMYA(J2-2:J2))117:             DUMMYC(J2-2:J2) = MATMUL(RMATI,DUMMYA(J2-2:J2))
121:          ENDDO118:          ENDDO
430: !      STOP427: !      STOP
431: 428: 
432:       END SUBROUTINE MINPERMDISTRBCOM429:       END SUBROUTINE MINPERMDISTRBCOM
433: 430: 
434: !     ----------------------------------------------------------------------------------------------431: !     ----------------------------------------------------------------------------------------------
435: 432: 
436:       SUBROUTINE ORIENTA(X, T1S, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2, NATOMS, Q2, DEBUG)433:       SUBROUTINE ORIENTA(X, T1S, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2, NATOMS, Q2, DEBUG)
437: 434: 
438: !     This subroutine puts the configuration, X, of an atomic  system into a standard alignment, T1.435: !     This subroutine puts the configuration, X, of an atomic  system into a standard alignment, T1.
439: !436: !
440:       USE KEY, ONLY: EFIELDT, MULTISITEPYT437:       USE KEY, ONLY: EFIELDT
441: 438: 
442:       IMPLICIT NONE439:       IMPLICIT NONE
443:       INTEGER          :: NATOMS, I, J, J1, J2, JMAX1, JMAX2, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2440:       INTEGER          :: NATOMS, I, J, J1, J2, JMAX1, JMAX2, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2
444:       DOUBLE PRECISION :: X(3*NATOMS), T1(3*NATOMS)441:       DOUBLE PRECISION :: X(3*NATOMS), T1(3*NATOMS)
445:       DOUBLE PRECISION :: XS(3*NATOMS), T1S(3*NATOMS), T2S(3*NATOMS), DIST(NATOMS)442:       DOUBLE PRECISION :: XS(3*NATOMS), T1S(3*NATOMS), T2S(3*NATOMS), DIST(NATOMS)
446:       DOUBLE PRECISION :: AX(3), P(3), Q2(4), ROTM(3,3), ROTMINV(3,3)443:       DOUBLE PRECISION :: AX(3), P(3), Q2(4), ROTM(3,3), ROTMINV(3,3)
447:       DOUBLE PRECISION :: THETA, THETAH, COST, SINT, COSTH, SINTH, ST, FCT444:       DOUBLE PRECISION :: THETA, THETAH, COST, SINT, COSTH, SINTH, ST, FCT
448:       DOUBLE PRECISION :: CMX, CMY, CMZ, DMAX, DUMMY, PROJ, DMAX2, CUTOFF1, DTEMP445:       DOUBLE PRECISION :: CMX, CMY, CMZ, DMAX, DUMMY, PROJ, DMAX2, CUTOFF1, DTEMP
449:       LOGICAL          :: DEBUG446:       LOGICAL          :: DEBUG
450: 447: 
459:          CMX = CMX + X(J-2)456:          CMX = CMX + X(J-2)
460:          CMY = CMY + X(J-1)457:          CMY = CMY + X(J-1)
461:          CMZ = CMZ + X(J)458:          CMZ = CMZ + X(J)
462:       ENDDO459:       ENDDO
463: 460: 
464:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS461:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS
465:       DO I = 1, NATOMS462:       DO I = 1, NATOMS
466:          J = 3*I463:          J = 3*I
467:          XS(J-2) = X(J-2) - CMX464:          XS(J-2) = X(J-2) - CMX
468:          XS(J-1) = X(J-1) - CMY465:          XS(J-1) = X(J-1) - CMY
469: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)466:          XS(J)   = X(J)   - CMZ
470:          IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)          XS(J)   = X(J)   - CMZ 
471:       ENDDO467:       ENDDO
472: 468: 
473:       DMAX    = -1.D0469:       DMAX    = -1.D0
474:       NORBIT1 = 1470:       NORBIT1 = 1
475: 471: 
476:       IF (EFIELDT) THEN472:       IF (EFIELDT) THEN
477:          T1S(:) = XS(:)473:          T1S(:) = XS(:)
478:          GOTO 100474:          GOTO 100
479:       ENDIF475:       ENDIF
480: 476: 
683: 679: 
684:       ENDIF680:       ENDIF
685: 681: 
686:       END SUBROUTINE ROTATMXZ682:       END SUBROUTINE ROTATMXZ
687: 683: 
688: !     ----------------------------------------------------------------------------------------------684: !     ----------------------------------------------------------------------------------------------
689: 685: 
690:       SUBROUTINE ROTATM(T, X, Q2, NATOMS)686:       SUBROUTINE ROTATM(T, X, Q2, NATOMS)
691: !     takes the set of coordinates T for NATOMS number of atoms and returns X after rotation via the 687: !     takes the set of coordinates T for NATOMS number of atoms and returns X after rotation via the 
692: !     quaternion Q2 about the origin 688: !     quaternion Q2 about the origin 
693:       USE KEY, ONLY : EFIELDT, MULTISITEPYT689: 
694:       IMPLICIT NONE690:       IMPLICIT NONE
695: 691: 
696:       INTEGER          :: I, J, NATOMS692:       INTEGER          :: I, J, NATOMS
697:       DOUBLE PRECISION :: T(3*NATOMS), X(3*NATOMS), T1(1:3), Q2(4), RM(3,3) 693:       DOUBLE PRECISION :: T(3*NATOMS), X(3*NATOMS), T1(1:3), Q2(4), RM(3,3) 
698:       DOUBLE PRECISION :: CMX, CMY, CMZ694:       DOUBLE PRECISION :: CMX, CMY, CMZ
699: !695: !
700: !     Move centre of mass to the origin.696: !     Move centre of mass to the origin.
701: !697: !
702:       CMX = 0.D0; CMY = 0.D0; CMZ = 0.D0698:       CMX = 0.D0; CMY = 0.D0; CMZ = 0.D0
703:       DO I = 1, NATOMS699:       DO I = 1, NATOMS
704:          J = 3*I700:          J = 3*I
705:          CMX = CMX + T(J-2)701:          CMX = CMX + T(J-2)
706:          CMY = CMY + T(J-1)702:          CMY = CMY + T(J-1)
707:          CMZ = CMZ + T(J)703:          CMZ = CMZ + T(J)
708:       ENDDO704:       ENDDO
709:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS705:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS
710:       DO I = 1, NATOMS706:       DO I = 1, NATOMS
711:          J      = 3*I707:          J      = 3*I
712:          X(J-2) = T(J-2) - CMX708:          X(J-2) = T(J-2) - CMX
713:          X(J-1) = T(J-1) - CMY709:          X(J-1) = T(J-1) - CMY
714: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)710:          X(J)   = T(J)   - CMZ
715:          IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)          X(J)   = T(J)   - CMZ 
716:       ENDDO711:       ENDDO
717: 712: 
718: !     extract the rotation matrix corresponding to rotation via Q713: !     extract the rotation matrix corresponding to rotation via Q
719: 714: 
720:       CALL QROTMAT(Q2,RM)715:       CALL QROTMAT(Q2,RM)
721: 716: 
722:       DO I = 1, NATOMS717:       DO I = 1, NATOMS
723: 718: 
724:          J        = 3*I719:          J        = 3*I
725:          T1(1:3)  = MATMUL(RM(:,:), X(J-2:J))720:          T1(1:3)  = MATMUL(RM(:,:), X(J-2:J))
728:       ENDDO723:       ENDDO
729: 724: 
730:       END SUBROUTINE ROTATM725:       END SUBROUTINE ROTATM
731: 726: 
732: !     ----------------------------------------------------------------------------------------------727: !     ----------------------------------------------------------------------------------------------
733: 728: 
734:       SUBROUTINE MINDISTA(RA,RB,NATOMS,DIST,Q2,DEBUG)729:       SUBROUTINE MINDISTA(RA,RB,NATOMS,DIST,Q2,DEBUG)
735: 730: 
736: !     Returns DIST as the actual distance, rather than the squared distance731: !     Returns DIST as the actual distance, rather than the squared distance
737: 732: 
738:       USE KEY, ONLY: EFIELDT, MULTISITEPYT733:       USE KEY, ONLY: EFIELDT
739: 734: 
740:       IMPLICIT NONE735:       IMPLICIT NONE
741: 736: 
742:       INTEGER          :: J1, J2, J3, J4, NATOMS, NSIZE, JMIN, INFO737:       INTEGER          :: J1, J2, J3, J4, NATOMS, NSIZE, JMIN, INFO
743:       DOUBLE PRECISION :: RA(3*NATOMS), RB(3*NATOMS), DIST, QMAT(4,4), TEMPA(9*NATOMS), XM, YM, ZM, XP, YP, ZP738:       DOUBLE PRECISION :: RA(3*NATOMS), RB(3*NATOMS), DIST, QMAT(4,4), TEMPA(9*NATOMS), XM, YM, ZM, XP, YP, ZP
744:       DOUBLE PRECISION :: DIAG(4), MINV, Q2(4), CMXA, CMYA, CMZA, CMXB, CMYB, CMZB739:       DOUBLE PRECISION :: DIAG(4), MINV, Q2(4), CMXA, CMYA, CMZA, CMXB, CMYB, CMZB
745:       DOUBLE PRECISION :: R(3), P(3), RM(3,3), THETA, FCT740:       DOUBLE PRECISION :: R(3), P(3), RM(3,3), THETA, FCT
746:       DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)741:       DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)
747:       DOUBLE PRECISION :: ENERGY, VNEW(3*NATOMS), RMS, DUMMY742:       DOUBLE PRECISION :: ENERGY, VNEW(3*NATOMS), RMS, DUMMY
748:       LOGICAL          :: BULKT, PRESERVET, DEBUG743:       LOGICAL          :: BULKT, PRESERVET, DEBUG
758:          J2 = 3*(J1-1)753:          J2 = 3*(J1-1)
759:          CMXA = CMXA + XA(J2+1)754:          CMXA = CMXA + XA(J2+1)
760:          CMYA = CMYA + XA(J2+2)755:          CMYA = CMYA + XA(J2+2)
761:          CMZA = CMZA + XA(J2+3)756:          CMZA = CMZA + XA(J2+3)
762:       ENDDO757:       ENDDO
763:       CMXA = CMXA/NSIZE; CMYA = CMYA/NSIZE; CMZA = CMZA/NSIZE758:       CMXA = CMXA/NSIZE; CMYA = CMYA/NSIZE; CMZA = CMZA/NSIZE
764:       DO J1 = 1, NSIZE759:       DO J1 = 1, NSIZE
765:          J2 = 3*(J1-1)760:          J2 = 3*(J1-1)
766:          XA(J2+1) = XA(J2+1) - CMXA761:          XA(J2+1) = XA(J2+1) - CMXA
767:          XA(J2+2) = XA(J2+2) - CMYA762:          XA(J2+2) = XA(J2+2) - CMYA
768: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)763:          XA(J2+3) = XA(J2+3) - CMZA
769:          IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)          XA(J2+3) = XA(J2+3) - CMZA 
770:       ENDDO764:       ENDDO
771: 765: 
772:       CMXB = 0.0D0; CMYB = 0.0D0; CMZB = 0.0D0766:       CMXB = 0.0D0; CMYB = 0.0D0; CMZB = 0.0D0
773:       DO J1 = 1, NSIZE767:       DO J1 = 1, NSIZE
774:          J2 = 3*(J1-1)768:          J2 = 3*(J1-1)
775:          CMXB = CMXB + XB(J2+1)769:          CMXB = CMXB + XB(J2+1)
776:          CMYB = CMYB + XB(J2+2)770:          CMYB = CMYB + XB(J2+2)
777:          CMZB = CMZB + XB(J2+3)771:          CMZB = CMZB + XB(J2+3)
778:       ENDDO772:       ENDDO
779:       CMXB = CMXB/NSIZE; CMYB = CMYB/NSIZE; CMZB = CMZB/NSIZE773:       CMXB = CMXB/NSIZE; CMYB = CMYB/NSIZE; CMZB = CMZB/NSIZE
780:       DO J1 = 1, NSIZE774:       DO J1 = 1, NSIZE
781:          J2 = 3*(J1-1)775:          J2 = 3*(J1-1)
782:          XB(J2+1) = XB(J2+1) - CMXB776:          XB(J2+1) = XB(J2+1) - CMXB
783:          XB(J2+2) = XB(J2+2) - CMYB777:          XB(J2+2) = XB(J2+2) - CMYB
784: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)778:          XB(J2+3) = XB(J2+3) - CMZB
785:          IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)          XB(J2+3) = XB(J2+3) - CMZB 
786:       ENDDO779:       ENDDO
787: 780: 
788:       QMAT(1:4,1:4) = 0.0D0781:       QMAT(1:4,1:4) = 0.0D0
789: 782: 
790:       DO J1 = 1, NSIZE783:       DO J1 = 1, NSIZE
791:          J2 = 3*(J1-1)784:          J2 = 3*(J1-1)
792:          XM = XA(J2+1) - XB(J2+1)785:          XM = XA(J2+1) - XB(J2+1)
793:          YM = XA(J2+2) - XB(J2+2)786:          YM = XA(J2+2) - XB(J2+2)
794:          ZM = XA(J2+3) - XB(J2+3)787:          ZM = XA(J2+3) - XB(J2+3)
795:          XP = XA(J2+1) + XB(J2+1)788:          XP = XA(J2+1) + XB(J2+1)
887:         Q2(1) = QMAT(1,JMIN); Q2(2) = QMAT(2,JMIN); Q2(3) = QMAT(3,JMIN); Q2(4) = QMAT(4,JMIN)880:         Q2(1) = QMAT(1,JMIN); Q2(2) = QMAT(2,JMIN); Q2(3) = QMAT(3,JMIN); Q2(4) = QMAT(4,JMIN)
888: 881: 
889:       ENDIF882:       ENDIF
890: 883: 
891:       DIST = SQRT(MINV)884:       DIST = SQRT(MINV)
892: 885: 
893:       DO J1 = 1, NATOMS886:       DO J1 = 1, NATOMS
894:          J2 = 3*(J1-1)887:          J2 = 3*(J1-1)
895:          RB(J2+1) = RB(J2+1) - CMXB888:          RB(J2+1) = RB(J2+1) - CMXB
896:          RB(J2+2) = RB(J2+2) - CMYB889:          RB(J2+2) = RB(J2+2) - CMYB
897: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)890:          RB(J2+3) = RB(J2+3) - CMZB
898:          IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)          RB(J2+3) = RB(J2+3) - CMZB 
899:       ENDDO891:       ENDDO
900: 892: 
901:       CALL NEWROTGEOMA(NATOMS,RB,Q2,RM,CMXA,CMYA,CMZA)893:       CALL NEWROTGEOMA(NATOMS,RB,Q2,RM,CMXA,CMYA,CMZA)
902: 894: 
903:       DEALLOCATE(XA,XB)895:       DEALLOCATE(XA,XB)
904: 896: 
905:       END SUBROUTINE MINDISTA897:       END SUBROUTINE MINDISTA
906: 898: 
907: !     ----------------------------------------------------------------------------------------------899: !     ----------------------------------------------------------------------------------------------
908: 900: 


r31852/multisitepy.f90 2017-01-30 12:30:11.080044779 +0000 r31851/multisitepy.f90 2017-01-30 12:30:15.428103911 +0000
 12: !     AAtoSites -  12: !     AAtoSites - 
 13: !     TAKESTEPMULTISITEPY -  13: !     TAKESTEPMULTISITEPY - 
 14:  14: 
 15: !-------------------------------------------------------! 15: !-------------------------------------------------------!
 16: ! This is the older version of the multisite subroutine ! 16: ! This is the older version of the multisite subroutine !
 17: !-------------------------------------------------------! 17: !-------------------------------------------------------!
 18:  18: 
 19:       SUBROUTINE MULTISITEPY (X, G, ENERGY, GTEST) 19:       SUBROUTINE MULTISITEPY (X, G, ENERGY, GTEST)
 20:  20: 
 21:       USE COMMONS, ONLY: NATOMS 21:       USE COMMONS, ONLY: NATOMS
 22:       USE KEY, ONLY :  NPYSITE, PYA1, PYA2, PYSIGNOT, PYEPSNOT, RADIFT, FROZEN, EFIELDT  22:       USE KEY, ONLY :  NPYSITE, PYA1, PYA2, PYSIGNOT, PYEPSNOT, RADIFT, FROZEN 
 23:       USE PYMODULE, ONLY : pi, twopi 23:       USE PYMODULE, ONLY : pi, twopi
 24:       IMPLICIT NONE 24:       IMPLICIT NONE
 25:  25: 
 26:       INTEGER          :: I, I1, J, J1, J2, J3, J4, J5, J6, K1, K2, REALNATOMS, OFFSET 26:       INTEGER          :: I, J, J1, J2, J3, J4, J5, J6, K1, K2, REALNATOMS, OFFSET
 27:       DOUBLE PRECISION :: X(3*NATOMS), G(3*NATOMS), energy_contrib2, grad_contrib2 27:       DOUBLE PRECISION :: X(3*NATOMS), G(3*NATOMS)
 28:       DOUBLE PRECISION :: PST(NATOMS,3), EZRI1(3,3), EZRI2(3,3), EZRJ1(3,3), EZRJ2(3,3) 28:       DOUBLE PRECISION :: PST(NATOMS,3), EZRI1(3,3), EZRI2(3,3), EZRJ1(3,3), EZRJ2(3,3)
 29:       DOUBLE PRECISION :: RI(3), RJ(3), RIJ(3), NR(3), RIJSQ, ABSRIJ, P(3), THETA, THETA2, CT, ST 29:       DOUBLE PRECISION :: RI(3), RJ(3), RIJ(3), NR(3), RIJSQ, ABSRIJ, P(3), THETA, THETA2, CT, ST
 30:       DOUBLE PRECISION :: I3(3,3), AE1(3,3), BE1(3,3), AE2(3,3), BE2(3,3), APB(3,3), APBINV(3,3) 30:       DOUBLE PRECISION :: I3(3,3), AE1(3,3), BE1(3,3), AE2(3,3), BE2(3,3), APB(3,3), APBINV(3,3)
 31:       DOUBLE PRECISION :: RSTI(3), RSTJ(3), FCNT1, FCNT2, SRTFI1, SRTFI2, FMIN, LAMDAC, ENERGY 31:       DOUBLE PRECISION :: RSTI(3), RSTJ(3), FCNT1, FCNT2, SRTFI1, SRTFI2, FMIN, LAMDAC, ENERGY
 32:       DOUBLE PRECISION :: RHO1, RHO1SQ, RHO16, RHO112, RHO2, RHO2SQ, RHO26 32:       DOUBLE PRECISION :: RHO1, RHO1SQ, RHO16, RHO112, RHO2, RHO2SQ, RHO26
 33:       DOUBLE PRECISION :: FCTR1, FCTR2, DV1DF1, DV2DF2, DV1DR, DV2DR 33:       DOUBLE PRECISION :: FCTR1, FCTR2, DV1DF1, DV2DF2, DV1DR, DV2DR
 34:       DOUBLE PRECISION :: DRDPI1, DRDPI2, DRDPI3, DRDPJ1, DRDPJ2, DRDPJ3  34:       DOUBLE PRECISION :: DRDPI1, DRDPI2, DRDPI3, DRDPJ1, DRDPJ2, DRDPJ3 
 35:       DOUBLE PRECISION :: DF1PI1, DF1PI2, DF1PI3, DF2PI1, DF2PI2, DF2PI3 35:       DOUBLE PRECISION :: DF1PI1, DF1PI2, DF1PI3, DF2PI1, DF2PI2, DF2PI3
 36:       DOUBLE PRECISION :: DF1PJ1, DF1PJ2, DF1PJ3, DF2PJ1, DF2PJ2, DF2PJ3  36:       DOUBLE PRECISION :: DF1PJ1, DF1PJ2, DF1PJ3, DF2PJ1, DF2PJ2, DF2PJ3 
 37:       DOUBLE PRECISION :: RMI(3,3), RMJ(3,3), E(3,3) 37:       DOUBLE PRECISION :: RMI(3,3), RMJ(3,3), E(3,3)
 55:       ENDDO 55:       ENDDO
 56:  56: 
 57:  57: 
 58:       CALL DEFMSPY(PST) 58:       CALL DEFMSPY(PST)
 59:  59: 
 60:       IF (GTEST) THEN 60:       IF (GTEST) THEN
 61:  61: 
 62:          ENERGY = 0.D0 62:          ENERGY = 0.D0
 63:          G(:)   = 0.D0 63:          G(:)   = 0.D0
 64:  64: 
 65:  
 66:          DO J1 = 1, REALNATOMS - 1 65:          DO J1 = 1, REALNATOMS - 1
 67:  66: 
 68:             J3      = 3*J1 67:             J3      = 3*J1
 69:             J5      = OFFSET + J3 68:             J5      = OFFSET + J3
 70:             RI      = X(J3-2:J3) 69:             RI      = X(J3-2:J3)
 71:             P       = X(J5-2:J5) 70:             P       = X(J5-2:J5)
 72:  71: 
 73:             angle=sqrt(dot_product(P,P)) 72:             angle=sqrt(dot_product(P,P))
 74:             if(angle>twopi) then 73:             if(angle>twopi) then
 75: ! normalise angle-axis coordinates 74: ! normalise angle-axis coordinates
303: 302: 
304:             ENDDO303:             ENDDO
305: 304: 
306: !     END OUTER LOOP OVER PARTICLES305: !     END OUTER LOOP OVER PARTICLES
307: 306: 
308:          ENDDO307:          ENDDO
309: 308: 
310:          ENERGY = 4.D0*PYEPSNOT*ENERGY309:          ENERGY = 4.D0*PYEPSNOT*ENERGY
311:          G      = 24.D0*PYEPSNOT*G310:          G      = 24.D0*PYEPSNOT*G
312: 311: 
313:     do i1 = 1, natoms/2 
314: !sf344> gravity field with constant gradient, and an r^-12 repulsion around z=0 
315:       if(efieldt) then 
316:         energy_contrib2=0.0 
317:         grad_contrib2=0.0 ! this is a scalar, we have only contribution along the z direction 
318:         call gravity(x(3*i1),energy_contrib2,grad_contrib2) 
319:         energy = energy + energy_contrib2 
320:         g(3*i1)=g(3*i1)+grad_contrib2 
321:       end if 
322:     end do 
323:  
324:       ELSE312:       ELSE
325: 313: 
326:          ENERGY = 0.D0314:          ENERGY = 0.D0
327: 315: 
328:          DO J1 = 1, REALNATOMS - 1316:          DO J1 = 1, REALNATOMS - 1
329: 317: 
330:             J3      = 3*J1318:             J3      = 3*J1
331:             J5      = OFFSET + J3319:             J5      = OFFSET + J3
332:             RI      = X(J3-2:J3)320:             RI      = X(J3-2:J3)
333:             P       = X(J5-2:J5)321:             P       = X(J5-2:J5)
466: 454: 
467:                ENDDO455:                ENDDO
468: 456: 
469: !     CALCULATE PY POTENTIAL ENERGY457: !     CALCULATE PY POTENTIAL ENERGY
470: 458: 
471:             ENDDO459:             ENDDO
472: 460: 
473:          ENDDO461:          ENDDO
474: 462: 
475:          ENERGY = 4.D0*PYEPSNOT*ENERGY463:          ENERGY = 4.D0*PYEPSNOT*ENERGY
476:     do i1 = 1, natoms/2 
477: !sf344> gravity field with constant gradient, and an r^-12 repulsion around z=0 
478:       if(efieldt) then 
479:         energy_contrib2=0.0 
480:         grad_contrib2=0.0 ! this is a scalar, we have only contribution along the z direction 
481:         call gravity(x(3*i1),energy_contrib2,grad_contrib2) 
482:         energy = energy + energy_contrib2 
483:         g(3*i1)=g(3*i1)+grad_contrib2 
484:       end if 
485:     end do 
486:  
487: 464: 
488:       ENDIF465:       ENDIF
489:       END SUBROUTINE MULTISITEPY 466:       END SUBROUTINE MULTISITEPY 
490: 467: 
491: 468: 
492: !---------------------------------!469: !---------------------------------!
493: ! DEFMSPY - Called by MULTISITEPY !470: ! DEFMSPY - Called by MULTISITEPY !
494: !---------------------------------!471: !---------------------------------!
495: ! This subroutine will define where the multiple PY sites are supposed to go.472: ! This subroutine will define where the multiple PY sites are supposed to go.
496: ! Because this is all hardcoded and is for MULTISITEPY, I think this is obsolete.473: ! Because this is all hardcoded and is for MULTISITEPY, I think this is obsolete.
634:       IF(.NOT.ALLOCATED(SMABF)) ALLOCATE(SMABF(NPYSITE,3,3))611:       IF(.NOT.ALLOCATED(SMABF)) ALLOCATE(SMABF(NPYSITE,3,3))
635:       IF(.NOT.ALLOCATED(SMRDET)) ALLOCATE(SMRDET(NPYSITE))612:       IF(.NOT.ALLOCATED(SMRDET)) ALLOCATE(SMRDET(NPYSITE))
636:       IF(.NOT.ALLOCATED(SMADET)) ALLOCATE(SMADET(NPYSITE))613:       IF(.NOT.ALLOCATED(SMADET)) ALLOCATE(SMADET(NPYSITE))
637:       IF(.NOT.ALLOCATED(ELLMAT)) ALLOCATE(ELLMAT(NPYSITE,3,3))614:       IF(.NOT.ALLOCATED(ELLMAT)) ALLOCATE(ELLMAT(NPYSITE,3,3))
638:       IF(.NOT.ALLOCATED(SITECOORDS)) ALLOCATE(SITECOORDS(NPYSITE,3))615:       IF(.NOT.ALLOCATED(SITECOORDS)) ALLOCATE(SITECOORDS(NPYSITE,3))
639:       IF(.NOT.ALLOCATED(XEND)) ALLOCATE(XEND(2,NPYSITE,3))616:       IF(.NOT.ALLOCATED(XEND)) ALLOCATE(XEND(2,NPYSITE,3))
640: 617: 
641:       IF(.NOT.ALLOCATED(RBSITE)) ALLOCATE(RBSITE(NRBSITES,3))618:       IF(.NOT.ALLOCATED(RBSITE)) ALLOCATE(RBSITE(NRBSITES,3))
642:       RADIFTARRAY(:) = .TRUE.619:       RADIFTARRAY(:) = .TRUE.
643:           vecsbf(1)=0.0D0620:           vecsbf(1)=0.0D0
644:           vecsbf(2)=0.0D0621:           vecsbf(2)=1.0D0
645:           vecsbf(3)=1.0D0622:           vecsbf(3)=0.0D0
646: !    IF(LJSITECOORDST) THEN623: !    IF(LJSITECOORDST) THEN
647: !        vecsbf(:)=LJSITECOORDS(:)624: !        vecsbf(:)=LJSITECOORDS(:)
648: !        WRITE(*,'(A,3F8.3)') ' repulsive LJ site coordinates will be ', LJSITECOORDS(:)625: !        WRITE(*,'(A,3F8.3)') ' repulsive LJ site coordinates will be ', LJSITECOORDS(:)
649: !    END IF626: !    END IF
650: 627: 
651:       ! Begin loop lines in input. Each line corresponds to a site in the molecule.628:       ! Begin loop lines in input. Each line corresponds to a site in the molecule.
652:       DO J1=1,NPYSITE629:       DO J1=1,NPYSITE
653:             ! Read line as per syntax above.630:             ! Read line as per syntax above.
654:             READ(299,*) LABEL, PST(J1,1), PST(J1,2), PST(J1,3), &631:             READ(299,*) LABEL, PST(J1,1), PST(J1,2), PST(J1,3), &
655:              & DUMMYLABEL1, pyrepp(J1), pyattrp(J1), &632:              & DUMMYLABEL1, pyrepp(J1), pyattrp(J1), &
689:             SMRDET(J1) = 1.D0 / (SMRBF(J1,1,1)*SMRBF(J1,2,2)*SMRBF(J1,3,3))666:             SMRDET(J1) = 1.D0 / (SMRBF(J1,1,1)*SMRBF(J1,2,2)*SMRBF(J1,3,3))
690:             SMADET(J1) = 1.D0 / (SMABF(J1,1,1)*SMABF(J1,2,2)*SMABF(J1,3,3))667:             SMADET(J1) = 1.D0 / (SMABF(J1,1,1)*SMABF(J1,2,2)*SMABF(J1,3,3))
691: 668: 
692:             ! Rewrite the shape matrices as they appear in the body frame669:             ! Rewrite the shape matrices as they appear in the body frame
693:             CALL RMDRVT(OST(J1,:), RMAT, DRMAT1, DRMAT2, DRMAT3, .FALSE.)670:             CALL RMDRVT(OST(J1,:), RMAT, DRMAT1, DRMAT2, DRMAT3, .FALSE.)
694:             SMRBF(J1,:,:) = MATMUL(RMAT,(MATMUL(SMRBF(J1,:,:),TRANSPOSE(RMAT))))671:             SMRBF(J1,:,:) = MATMUL(RMAT,(MATMUL(SMRBF(J1,:,:),TRANSPOSE(RMAT))))
695:             SMABF(J1,:,:) = MATMUL(RMAT,(MATMUL(SMABF(J1,:,:),TRANSPOSE(RMAT))))672:             SMABF(J1,:,:) = MATMUL(RMAT,(MATMUL(SMABF(J1,:,:),TRANSPOSE(RMAT))))
696: 673: 
697: ! generate coordinates for the end of the largest ellipsoid axis674: ! generate coordinates for the end of the largest ellipsoid axis
698:         ! sf344> reset RBSITES to the centre of mass for each ellipsoid. Forget about RBSITE for now.675:         ! sf344> reset RBSITES to the centre of mass for each ellipsoid. Forget about RBSITE for now.
699: ! sf344> try to define sites for axially symmetric particles...676:             xend(1,J1,:)=PST(J1,:)!+ELLST1(J1,2)*MATMUL(RMAT,vecsbf)    ! vecsbf: (1,0,0) in the body frame of ellipsoid
700:             xend(1,J1,:)=PST(J1,:)+MATMUL(RMAT,vecsbf) !+ELLST1(J1,2)*MATMUL(RMAT,vecsbf)    ! vecsbf: (1,0,0) in the body frame of ellipsoid677:             xend(2,J1,:)=PST(J1,:)!-ELLST1(J1,2)*MATMUL(RMAT,vecsbf)    ! vecsbf: (1,0,0) in the body frame of ellipsoid
701:             xend(2,J1,:)=PST(J1,:)-MATMUL(RMAT,vecsbf) !-ELLST1(J1,2)*MATMUL(RMAT,vecsbf)    ! vecsbf: (1,0,0) in the body frame of ellipsoid 
702: 678: 
703:       END DO ! loop over sites679:       END DO ! loop over sites
704: 680: 
705:       J1=0681:       J1=0
706:       J2=0682:       J2=0
707:       DO 683:       DO 
708:         J1=J1+2684:         J1=J1+2
709:         J2=J2+1685:         J2=J2+1
710:          RBSITE(J1-1,:)=XEND(1,J2,:)686:          RBSITE(J1-1,:)=XEND(1,J2,:)
711:          RBSITE(J1,:)=XEND(2,J2,:)687:          RBSITE(J1,:)=XEND(2,J2,:)
860:       ! LJGSITEEPS, LJGSITESIGMA = cf DEFINELJMULTISITES above836:       ! LJGSITEEPS, LJGSITESIGMA = cf DEFINELJMULTISITES above
861:       ! [CUTOFF]837:       ! [CUTOFF]
862:       ! PARAMONOVPBCX = logical if there is a PBC in the x direction838:       ! PARAMONOVPBCX = logical if there is a PBC in the x direction
863:       ! BOXLX = size of PBC box in x direction in absolute units839:       ! BOXLX = size of PBC box in x direction in absolute units
864:       ! PARAMONOVCUTOFF = logical if cutoff is used840:       ! PARAMONOVCUTOFF = logical if cutoff is used
865:       ! PCUTOFF = cutoff distance in absolute units841:       ! PCUTOFF = cutoff distance in absolute units
866:       USE COMMONS, ONLY: NATOMS 842:       USE COMMONS, ONLY: NATOMS 
867:       USE KEY, ONLY: NPYSITE, PYSIGNOT, PYEPSNOT, FROZEN, &843:       USE KEY, ONLY: NPYSITE, PYSIGNOT, PYEPSNOT, FROZEN, &
868:             & LJGSITET, LJGSITEEPS, LJGSITESIGMA, &844:             & LJGSITET, LJGSITEEPS, LJGSITESIGMA, &
869:             & PARAMONOVPBCX, PARAMONOVPBCY, PARAMONOVPBCZ, &845:             & PARAMONOVPBCX, PARAMONOVPBCY, PARAMONOVPBCZ, &
870:             & PARAMONOVCUTOFF, PCUTOFF, EFIELDT, BULKT846:             & PARAMONOVCUTOFF, PCUTOFF
871:       USE PYMODULE, ONLY : BOXLX, BOXLY, BOXLZ847:       USE PYMODULE, ONLY : BOXLX, BOXLY, BOXLZ
872: 848: 
873:       ! PST, ..., SMABF = cf DEFINEPYMULTISITES above849:       ! PST, ..., SMABF = cf DEFINEPYMULTISITES above
874:       ! NLJSITE, ..., ljattrp = cf DEFINELJMULTISITES above850:       ! NLJSITE, ..., ljattrp = cf DEFINELJMULTISITES above
875:       ! LONGRSAX, LONGASAX = cf DEFINEPYMULTISITES above851:       ! LONGRSAX, LONGASAX = cf DEFINEPYMULTISITES above
876:       USE PYMODULE, ONLY : PST, OST, RADIFTARRAY, pi, twopi, &852:       USE PYMODULE, ONLY : PST, OST, RADIFTARRAY, pi, twopi, &
877:             & pyrepp, pyattrp, LONGRSAX, LONGASAX, SMRBF, SMABF, &853:             & pyrepp, pyattrp, LONGRSAX, LONGASAX, SMRBF, SMABF, &
878:             & NLJSITE, LJGSITECOORDS, ljrepp, ljattrp, &854:             & NLJSITE, LJGSITECOORDS, ljrepp, ljattrp, &
879:             & SMRDET, SMADET855:             & SMRDET, SMADET
880:  856:  
993:       ! K = index for looping over images969:       ! K = index for looping over images
994:       ! KOUT = marker for closest image970:       ! KOUT = marker for closest image
995:       ! PBCCOUNT = number of images to test for closeness971:       ! PBCCOUNT = number of images to test for closeness
996:       ! ROUT = holder for putative min for g_1 (or max of RHO1,2)972:       ! ROUT = holder for putative min for g_1 (or max of RHO1,2)
997:       ! RJST[image][xyz] = "r_ij" sotre holding image locations973:       ! RJST[image][xyz] = "r_ij" sotre holding image locations
998:       ! LOUT = lambda out for lambda_c val from putative g_1 min974:       ! LOUT = lambda out for lambda_c val from putative g_1 min
999:       LOGICAL :: PBC975:       LOGICAL :: PBC
1000:       INTEGER :: K, KOUT, PBCCOUNT976:       INTEGER :: K, KOUT, PBCCOUNT
1001:       DOUBLE PRECISION :: ROUT, RJST(8,3), LOUT977:       DOUBLE PRECISION :: ROUT, RJST(8,3), LOUT
1002: 978: 
1003:       INTEGER :: I1 
1004:       DOUBLE PRECISION :: energy_contrib2,grad_contrib2 
1005: 979: 
1006:       !-------------------------!980:       !-------------------------!
1007:       ! Initialize variables #B !981:       ! Initialize variables #B !
1008:       !-------------------------!982:       !-------------------------!
1009: 983: 
1010:       ! NATOMS is the number of lines in the coords file. Because one line specifies984:       ! NATOMS is the number of lines in the coords file. Because one line specifies
1011:       ! the mol's position and the next specifies the components of p for a molecule,985:       ! the mol's position and the next specifies the components of p for a molecule,
1012:       ! the real number of atoms is NATOMS/2. Set up the offset in indices.986:       ! the real number of atoms is NATOMS/2. Set up the offset in indices.
1013:       ! Initialize PBC check.987:       ! Initialize PBC check.
1014:       REALNATOMS = NATOMS/2988:       REALNATOMS = NATOMS/2
1040:             GCUT = PYSIGNOT/(PCUTOFF+PYSIGNOT)1014:             GCUT = PYSIGNOT/(PCUTOFF+PYSIGNOT)
1041:             GCUT2 = GCUT*GCUT1015:             GCUT2 = GCUT*GCUT
1042:             GCUT6 = GCUT2**31016:             GCUT6 = GCUT2**3
1043:             GCUT12 = GCUT6**21017:             GCUT12 = GCUT6**2
1044:       ENDIF ! [CUTOFF]1018:       ENDIF ! [CUTOFF]
1045: 1019: 
1046:       ! Initialize energies (and gradient if required)1020:       ! Initialize energies (and gradient if required)
1047:       ENERGY = 0.D01021:       ENERGY = 0.D0
1048:       IF(GTEST) G(:)   = 0.D01022:       IF(GTEST) G(:)   = 0.D0
1049: 1023: 
1050:  
1051:       ! [LJGSITE]1024:       ! [LJGSITE]
1052:       IF(LJGSITET) THEN1025:       IF(LJGSITET) THEN
1053:           LJENERGY = 0.D01026:           LJENERGY = 0.D0
1054:           IF(GTEST) LJG(:) = 0.D01027:           IF(GTEST) LJG(:) = 0.D0
1055:       ENDIF ! [LJGSITE]1028:       ENDIF ! [LJGSITE]
1056: 1029: 
1057:     ! There are a number of calculations that depend only on the positions, orientations,1030:     ! There are a number of calculations that depend only on the positions, orientations,
1058:     ! and parameters of each site, so it is more efficient to do all those once at the 1031:     ! and parameters of each site, so it is more efficient to do all those once at the 
1059:     ! beginning rather than repeat them inside the site loops.1032:     ! beginning rather than repeat them inside the site loops.
1060:     ! Loop over all mols1033:     ! Loop over all mols
1765:       ENDDO ! outer loop over molecules1738:       ENDDO ! outer loop over molecules
1766: 1739: 
1767:       !-------------!1740:       !-------------!
1768:       ! Finalize #N !1741:       ! Finalize #N !
1769:       !-------------!1742:       !-------------!
1770:       1743:       
1771:       ! Tack on the final factors to the energy and gradient as per eq 24 & 301744:       ! Tack on the final factors to the energy and gradient as per eq 24 & 30
1772:       ENERGY = 4.D0*PYEPSNOT*ENERGY1745:       ENERGY = 4.D0*PYEPSNOT*ENERGY
1773:       IF(GTEST) G = 24.D0*PYEPSNOT*G1746:       IF(GTEST) G = 24.D0*PYEPSNOT*G
1774: 1747: 
1775:      ![GRAVITY] 
1776:      do i1=1,realnatoms 
1777: !sf344> gravity field with constant gradient, and an r^-12 repulsion around z=-10 (for backwards compatibility with existing databases) 
1778:       if(efieldt) then 
1779:         energy_contrib2=0.0 
1780:         grad_contrib2=0.0 ! this is a scalar, we have only contribution along the z direction 
1781:         call gravity(x(3*i1),energy_contrib2,grad_contrib2) 
1782:         energy = energy + energy_contrib2 
1783:         g(3*i1)=g(3*i1)+grad_contrib2 
1784:       end if 
1785:      end do 
1786:  
1787:       ! [LJGSITE]1748:       ! [LJGSITE]
1788:       ! Add the LJ contributions to the PY contributions1749:       ! Add the LJ contributions to the PY contributions
1789:       IF(LJGSITET) THEN1750:       IF(LJGSITET) THEN
1790:             ! Tack on final factors to energy and gradient1751:             ! Tack on final factors to energy and gradient
1791:             LJENERGY = 4.D0*LJGSITEEPS*LJENERGY1752:             LJENERGY = 4.D0*LJGSITEEPS*LJENERGY
1792:             IF(GTEST) LJG = 4.D0*LJGSITEEPS*LJG 1753:             IF(GTEST) LJG = 4.D0*LJGSITEEPS*LJG 
1793: 1754: 
1794:             ENERGY = ENERGY + LJENERGY1755:             ENERGY = ENERGY + LJENERGY
1795:             IF(GTEST) G = G + LJG1756:             IF(GTEST) G = G + LJG
1796:       ENDIF ! [LJGSITE]1757:       ENDIF ! [LJGSITE]
2325: COFA(2,1) = COFA(1,2)2286: COFA(2,1) = COFA(1,2)
2326: COFA(2,2) = A(3,3)*A(1,1) - A(1,3)*A(3,1)2287: COFA(2,2) = A(3,3)*A(1,1) - A(1,3)*A(3,1)
2327: COFA(2,3) = A(3,1)*A(1,2) - A(1,1)*A(3,2)2288: COFA(2,3) = A(3,1)*A(1,2) - A(1,1)*A(3,2)
2328: COFA(3,1) = COFA(1,3)2289: COFA(3,1) = COFA(1,3)
2329: COFA(3,2) = COFA(2,3)2290: COFA(3,2) = COFA(2,3)
2330: COFA(3,3) = A(1,1)*A(2,2) - A(2,1)*A(1,2)2291: COFA(3,3) = A(1,1)*A(2,2) - A(2,1)*A(1,2)
2331: 2292: 
2332: RETURN2293: RETURN
2333: 2294: 
2334: END SUBROUTINE COFACTORS2295: END SUBROUTINE COFACTORS
2335:  
2336:  
2337: SUBROUTINE INERTIAPY(RMI, KBLOCK, TMASS) 
2338: ! sf344> first implementation of the inertia routine for multisite PY particles 
2339: !        applicable only for those centred around the origin 
2340: !        copied over from INERTIAPAP 
2341:  
2342:         IMPLICIT NONE 
2343:  
2344:         DOUBLE PRECISION, INTENT(IN)  :: RMI(3,3) 
2345:         DOUBLE PRECISION, INTENT(OUT) :: TMASS, KBLOCK(3,3) 
2346: ! LFCT: half the equilibrium pair distance squared 
2347: ! MOMENT: scalar moment of inertia of a sphere about any axis 
2348:         DOUBLE PRECISION :: LFCT, MOMENT 
2349:  
2350: ! Set the total mass as unity 
2351:         TMASS  = 1.D0 
2352: ! Set the value of LFCT to half the equilibrium pair distance squared 
2353:         LFCT = 1.44D0 
2354: ! Set the value of MOMENT using the formula I= 0.2*M*R*R 
2355:         MOMENT = 0.2D0*TMASS*LFCT 
2356:  
2357: ! Set the moment of inertia tensor to zero. 
2358: ! Off diagonal elements will remain zero 
2359:         KBLOCK(:,:) = 0.D0 
2360:  
2361: ! Set the diagonal elements to MOMENT 
2362:         KBLOCK(1,1) = MOMENT 
2363:         KBLOCK(2,2) = MOMENT 
2364:         KBLOCK(3,3) = MOMENT 
2365:  
2366:       END SUBROUTINE INERTIAPY 
2367:  


r31852/ncutils.f90 2017-01-30 12:30:08.772013392 +0000 r31851/ncutils.f90 2017-01-30 12:30:12.972070507 +0000
1300:           IMPLICIT NONE1300:           IMPLICIT NONE
1301: 1301: 
1302:           CHARACTER(LEN=132),INTENT(IN)         :: C1302:           CHARACTER(LEN=132),INTENT(IN)         :: C
1303:           CHARACTER(LEN=5),POINTER,DIMENSION(:) :: S1303:           CHARACTER(LEN=5),POINTER,DIMENSION(:) :: S
1304:           DOUBLE PRECISION,POINTER,DIMENSION(:) :: Q1304:           DOUBLE PRECISION,POINTER,DIMENSION(:) :: Q
1305:           DOUBLE PRECISION SITES(3*NTSITES), P(3), RM(3,3)1305:           DOUBLE PRECISION SITES(3*NTSITES), P(3), RM(3,3)
1306:           CHARACTER(LEN=2) ZSTRING(NATOMS)1306:           CHARACTER(LEN=2) ZSTRING(NATOMS)
1307: 1307: 
1308:           INTEGER :: J1308:           INTEGER :: J
1309:           DOUBLE PRECISION :: TMPCOORDS(9*NATOMS/2)1309:           DOUBLE PRECISION :: TMPCOORDS(9*NATOMS/2)
 1310: 
1310:           IF (STOCKT) THEN1311:           IF (STOCKT) THEN
1311:              WRITE(50,'(I6)') (NATOMS/2)1312:              WRITE(50,'(I6)') (NATOMS/2)
1312:              WRITE(50,'(1X,A)') TRIM(ADJUSTL(C))1313:              WRITE(50,'(1X,A)') TRIM(ADJUSTL(C))
1313:              DO J=1,(NATOMS/2)1314:              DO J=1,(NATOMS/2)
1314:                 WRITE(50,'(A5,1X,6F20.10)') S(J),Q(3*(J-1)+1),Q(3*(J-1)+2),Q(3*(J-1)+3), &1315:                 WRITE(50,'(A5,1X,6F20.10)') S(J),Q(3*(J-1)+1),Q(3*(J-1)+2),Q(3*(J-1)+3), &
1315:   &                             Q(3*((NATOMS/2)+J-1)+1),Q(3*((NATOMS/2)+J-1)+2),Q(3*((NATOMS/2)+J-1)+3)1316:   &                             Q(3*((NATOMS/2)+J-1)+1),Q(3*((NATOMS/2)+J-1)+2),Q(3*((NATOMS/2)+J-1)+3)
1316:              ENDDO1317:              ENDDO
1317: !          ELSE IF (RBAAT) THEN1318: !          ELSE IF (RBAAT) THEN
1318: !                WRITE(50,'(I6)') (NATOMS/2)1319: !                WRITE(50,'(I6)') (NATOMS/2)
1319: !                WRITE(50,'(1X,A)') TRIM(ADJUSTL(C))1320: !                WRITE(50,'(1X,A)') TRIM(ADJUSTL(C))


r31852/nnutils.f90 2017-01-30 12:30:09.000016492 +0000 r31851/nnutils.f90 2017-01-30 12:30:13.200073608 +0000
110:           DEVIATION = ABS( 100*( (NIMAGE+1)*DVEC/SEPARATION -1 ) )110:           DEVIATION = ABS( 100*( (NIMAGE+1)*DVEC/SEPARATION -1 ) )
111:           AVDEV     = SUM(DEVIATION)/(NIMAGE+1)111:           AVDEV     = SUM(DEVIATION)/(NIMAGE+1)
112:      END SUBROUTINE INTERNALIMAGEDISTRIBUTION112:      END SUBROUTINE INTERNALIMAGEDISTRIBUTION
113: 113: 
114: !--------------------------------------------------------------------------------------------------------114: !--------------------------------------------------------------------------------------------------------
115: 115: 
116:      SUBROUTINE MAKEIMAGE(EINITIAL,EFINAL,QQ,FINFIN) ! INTERPOLATES THE BAND116:      SUBROUTINE MAKEIMAGE(EINITIAL,EFINAL,QQ,FINFIN) ! INTERPOLATES THE BAND
117:           USE NEBDATA117:           USE NEBDATA
118:           USE SPFUNCTS118:           USE SPFUNCTS
119: !          USE KEYNEB,ONLY: NIMAGE119: !          USE KEYNEB,ONLY: NIMAGE
120:           USE KEY,ONLY: MORPHT, MSTEPS, GREATCIRCLET, GCIMAGE, GCCONV, GCUPDATE, GCSTEPS, INTEPSILON, MULTISITEPYT, PYGPERIODICT, PYT, &120:           USE KEY,ONLY: MORPHT, MSTEPS, GREATCIRCLET, GCIMAGE, GCCONV, GCUPDATE, GCSTEPS, INTEPSILON, &
121:                         BULKT, RBAAT, NTSITES, AMBERT, LOCALPERMDIST, NRBGROUP, RBGROUP, RBNINGROUP, NRBTRIES, NABT,TWOD, &121:                         BULKT, RBAAT, NTSITES, AMBERT, LOCALPERMDIST, NRBGROUP, RBGROUP, RBNINGROUP, NRBTRIES, NABT,TWOD, &
122:                         RIGIDBODY, AVOID_COLLISIONS, COLL_TOL, BULK_BOXVEC122:                         RIGIDBODY, AVOID_COLLISIONS, COLL_TOL, BULK_BOXVEC
123:           USE INTCOMMONS, ONLY : NNZ, KD, NINTC, DESMINT, DIHINFO, PREVDIH, BACKTCUTOFF, INTERPBACKTCUT, MINBACKTCUT, &123:           USE INTCOMMONS, ONLY : NNZ, KD, NINTC, DESMINT, DIHINFO, PREVDIH, BACKTCUTOFF, INTERPBACKTCUT, MINBACKTCUT, &
124:                                  CHICDNEB124:                                  CHICDNEB
125:           USE COMMONS, ONLY : ZSYM, NRBSITES, DEBUG125:           USE COMMONS, ONLY : ZSYM, NRBSITES, DEBUG
126:           USE MODCHARMM, ONLY : CHRMMT, ICINTERPT126:           USE MODCHARMM, ONLY : CHRMMT, ICINTERPT
127:           USE MODAMBER9, ONLY: AMBERICT, AMBICDNEBT, NICTOT127:           USE MODAMBER9, ONLY: AMBERICT, AMBICDNEBT, NICTOT
128:           USE PORFUNCS 128:           USE PORFUNCS 
129:           USE KEYNEB,ONLY: NIMAGE,XYZFILE,RBXYZFILE129:           USE KEYNEB,ONLY: NIMAGE,XYZFILE,RBXYZFILE
130:           USE KEY, ONLY: FILTH,FILTHSTR, MULTISITEPYT, ALIGNRBST, N_TO_ALIGN, TO_ALIGN, SLERPT 130:           USE KEY, ONLY: FILTH,FILTHSTR, MULTISITEPYT, ALIGNRBST, N_TO_ALIGN, TO_ALIGN, SLERPT 
369: 369: 
370: 370: 
371: 371: 
372: ! sf344> now we have the interpolated coordinates. For MULTISITEPY, check for overlap372: ! sf344> now we have the interpolated coordinates. For MULTISITEPY, check for overlap
373: !        and if there's any (very likely), move the offending bodies randomly around. 373: !        and if there's any (very likely), move the offending bodies randomly around. 
374:             IF(MULTISITEPYT) THEN    ! sn402 : added this IF statement374:             IF(MULTISITEPYT) THEN    ! sn402 : added this IF statement
375:                  DO I = 2, NIMAGE+1375:                  DO I = 2, NIMAGE+1
376:                     J = NOPT*(I-1)376:                     J = NOPT*(I-1)
377:                     CALL TAKESTEPMULTISITEPY(XYZ(J+1:J+NOPT))377:                     CALL TAKESTEPMULTISITEPY(XYZ(J+1:J+NOPT))
378:                  ENDDO378:                  ENDDO
379:            ELSE IF(PYT) THEN379:             ENDIF
380:               DO I = 2, NIMAGE+1 
381:                 J = NOPT*(I-1) 
382:                 CALL py_takestep(XYZ(J+1:J+NOPT)) 
383:               ENDDO 
384:            END IF 
385: 380: 
386:             DEALLOCATE(DELTAX)381:             DEALLOCATE(DELTAX)
387: 382: 
388:           ELSE   ! End of "else if(RBAAT .OR. RIGIDMOLECULEST)"383:           ELSE   ! End of "else if(RBAAT .OR. RIGIDMOLECULEST)"
389: 384: 
390:              ALLOCATE(DELTAX(NOPT))385:              ALLOCATE(DELTAX(NOPT))
391: 386: 
392:              IF (BULKT) THEN387:              IF (BULKT) THEN
393:                 DO K=1,NATOMS388:                 DO K=1,NATOMS
394:                    DELTAX(3*(K-1)+1)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1) &389:                    DELTAX(3*(K-1)+1)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1) &


r31852/OPTIM.F 2017-01-30 12:30:09.224019540 +0000 r31851/OPTIM.F 2017-01-30 12:30:13.428076709 +0000
391:          WRITE(*,'(A,G20.10)') ' OPTIM> Using z rotational ev shift=',SHIFTV391:          WRITE(*,'(A,G20.10)') ' OPTIM> Using z rotational ev shift=',SHIFTV
392:          SHIFTL(6)=SHIFTV392:          SHIFTL(6)=SHIFTV
393:       ELSE IF (TWISTT) THEN393:       ELSE IF (TWISTT) THEN
394:          WRITE(*,'(A,G20.10)') ' OPTIM> Using x,y,z trans and z rotation ev shift=',SHIFTV394:          WRITE(*,'(A,G20.10)') ' OPTIM> Using x,y,z trans and z rotation ev shift=',SHIFTV
395:          SHIFTL(1)=SHIFTV395:          SHIFTL(1)=SHIFTV
396:          SHIFTL(2)=SHIFTV396:          SHIFTL(2)=SHIFTV
397:          SHIFTL(3)=SHIFTV397:          SHIFTL(3)=SHIFTV
398:          SHIFTL(6)=SHIFTV398:          SHIFTL(6)=SHIFTV
399:          WRITE(*,'(A)') ' OPTIM> x,y translational and z rotational ev shifting'399:          WRITE(*,'(A)') ' OPTIM> x,y translational and z rotational ev shifting'
400:       ELSE IF (PULLT.OR.EFIELDT) THEN400:       ELSE IF (PULLT.OR.EFIELDT) THEN
401:          WRITE(*,'(A,G20.10)') ' OPTIM> Using x,y translational and z rotational ev shift=',SHIFTV401:          WRITE(*,'(A,G20.10)') ' OPTIM> Using x,y,z trans and z rotation ev shift=',SHIFTV
402:          SHIFTL(1)=SHIFTV402:          SHIFTL(1)=SHIFTV
403:          SHIFTL(2)=SHIFTV403:          SHIFTL(2)=SHIFTV
404: !         SHIFTL(3)=SHIFTV404:          SHIFTL(3)=SHIFTV
405:          SHIFTL(6)=SHIFTV405:          SHIFTL(6)=SHIFTV
406: !         WRITE(*,'(A)') ' OPTIM> x,y translational and z rotational ev shifting'406:          WRITE(*,'(A)') ' OPTIM> x,y translational and z rotational ev shifting'
407:       ELSE IF (RTEST) THEN407:       ELSE IF (RTEST) THEN
408:          IF (JZ.NE.0.0D0) THEN408:          IF (JZ.NE.0.0D0) THEN
409:             WRITE(*,'(A,G20.10)') ' OPTIM> Using trans/z rotation ev shift=',SHIFTV409:             WRITE(*,'(A,G20.10)') ' OPTIM> Using trans/z rotation ev shift=',SHIFTV
410:             SHIFTL(3)=SHIFTV410:             SHIFTL(3)=SHIFTV
411:             SHIFTL(6)=SHIFTV411:             SHIFTL(6)=SHIFTV
412:          ELSE412:          ELSE
413:             WRITE(*,'(A,G20.10)') ' OPTIM> Using z rotation ev shift=',SHIFTV413:             WRITE(*,'(A,G20.10)') ' OPTIM> Using z rotation ev shift=',SHIFTV
414:             SHIFTL(6)=SHIFTV414:             SHIFTL(6)=SHIFTV
415:          ENDIF415:          ENDIF
416: C     ELSE IF (NOSHIFT) THEN416: C     ELSE IF (NOSHIFT) THEN


r31852/path.f 2017-01-30 12:30:11.824055463 +0000 r31851/path.f 2017-01-30 12:30:15.664107115 +0000
809:                   CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),809:                   CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),
810:      1                          Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),810:      1                          Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),
811:      2                          CAPSCOORDS2,RAD,HEIGHT)811:      2                          CAPSCOORDS2,RAD,HEIGHT)
812:                   CALL CAPSIDIO(Q1(3*(J2-1)+1),Q1(3*(J2-1)+2),Q1(3*(J2-1)+3),812:                   CALL CAPSIDIO(Q1(3*(J2-1)+1),Q1(3*(J2-1)+2),Q1(3*(J2-1)+3),
813:      1                          Q1(3*(NATOMS/2+J2-1)+1),Q1(3*(NATOMS/2+J2-1)+2),Q1(3*(NATOMS/2+J2-1)+3),813:      1                          Q1(3*(NATOMS/2+J2-1)+1),Q1(3*(NATOMS/2+J2-1)+2),Q1(3*(NATOMS/2+J2-1)+3),
814:      2                          CAPSCOORDS1,RAD,HEIGHT)814:      2                          CAPSCOORDS1,RAD,HEIGHT)
815:                   DO J3=1,18815:                   DO J3=1,18
816:                       TEMP=TEMP+(CAPSCOORDS1(J3)-CAPSCOORDS2(J3))**2816:                       TEMP=TEMP+(CAPSCOORDS1(J3)-CAPSCOORDS2(J3))**2
817:                   ENDDO817:                   ENDDO
818:                ENDDO818:                ENDDO
819:             ELSE IF((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT) THEN819:             ELSE IF((PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT) THEN
820: C           calculate path lengths wrt. Cartesian coordinate of the centre and two coordinates along the symmetry axis only820: C           calculate path lengths wrt. Cartesian coordinate of the centre and two coordinates along the symmetry axis only
821:                CALL UNIAXGETPATHLENGTH(Q1,Q2,TEMP)821:                CALL UNIAXGETPATHLENGTH(Q1,Q2,TEMP)
822:             ELSE822:             ELSE
823:                DO J2=1,NOPT823:                DO J2=1,NOPT
824:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2824:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2
825:                ENDDO825:                ENDDO
826:             ENDIF826:             ENDIF
827:          ENDIF827:          ENDIF
828:          ! Compute the total pathlength from the TS to this frame828:          ! Compute the total pathlength from the TS to this frame
829:          PATHLENGTH(NSTEPPLUS+1-J1)=PATHLENGTH(NSTEPPLUS+2-J1)-SQRT(TEMP)829:          PATHLENGTH(NSTEPPLUS+1-J1)=PATHLENGTH(NSTEPPLUS+2-J1)-SQRT(TEMP)
932:                   CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),932:                   CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),
933:      1                          Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),933:      1                          Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),
934:      2                          CAPSCOORDS2,RAD,HEIGHT)934:      2                          CAPSCOORDS2,RAD,HEIGHT)
935:                   CALL CAPSIDIO(Q1(3*(J2-1)+1),Q1(3*(J2-1)+2),Q1(3*(J2-1)+3),935:                   CALL CAPSIDIO(Q1(3*(J2-1)+1),Q1(3*(J2-1)+2),Q1(3*(J2-1)+3),
936:      1                          Q1(3*(NATOMS/2+J2-1)+1),Q1(3*(NATOMS/2+J2-1)+2),Q1(3*(NATOMS/2+J2-1)+3),936:      1                          Q1(3*(NATOMS/2+J2-1)+1),Q1(3*(NATOMS/2+J2-1)+2),Q1(3*(NATOMS/2+J2-1)+3),
937:      2                          CAPSCOORDS1,RAD,HEIGHT)937:      2                          CAPSCOORDS1,RAD,HEIGHT)
938:                   DO J3=1,18938:                   DO J3=1,18
939:                       TEMP=TEMP+(CAPSCOORDS1(J3)-CAPSCOORDS2(J3))**2939:                       TEMP=TEMP+(CAPSCOORDS1(J3)-CAPSCOORDS2(J3))**2
940:                   ENDDO940:                   ENDDO
941:                ENDDO941:                ENDDO
942:             ELSE IF((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT) THEN942:             ELSE IF((PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT) THEN
943: C           calculate path lengths wrt. Cartesian coordinate of the centre and two coordinates along the symmetry axis only943: C           calculate path lengths wrt. Cartesian coordinate of the centre and two coordinates along the symmetry axis only
944:                CALL UNIAXGETPATHLENGTH(Q1,Q2,TEMP)944:                CALL UNIAXGETPATHLENGTH(Q1,Q2,TEMP)
945:             ELSE945:             ELSE
946:                DO J2=1,NOPT946:                DO J2=1,NOPT
947:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2947:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2
948:                ENDDO948:                ENDDO
949:             ENDIF949:             ENDIF
950:          ENDIF950:          ENDIF
951:          PATHLENGTH(NSTEPPLUS+1+J1)=PATHLENGTH(NSTEPPLUS+J1)+SQRT(TEMP)951:          PATHLENGTH(NSTEPPLUS+1+J1)=PATHLENGTH(NSTEPPLUS+J1)+SQRT(TEMP)
952:          IF (ZSYMSAVE(1:1).EQ.'W') THEN952:          IF (ZSYMSAVE(1:1).EQ.'W') THEN


r31852/potential.f 2017-01-30 12:30:12.052057995 +0000 r31851/potential.f 2017-01-30 12:30:15.892110215 +0000
1830:             ! WRITE(*,*) VNEW(J1)1830:             ! WRITE(*,*) VNEW(J1)
1831:             ! END DO1831:             ! END DO
1832:             ! STOP1832:             ! STOP
1833:             IF (SSTEST) THEN1833:             IF (SSTEST) THEN
1834:                CALL PYGPERIODICSECDER(COORDS,SSTEST)1834:                CALL PYGPERIODICSECDER(COORDS,SSTEST)
1835:             END IF1835:             END IF
1836:             IF (PTEST) THEN1836:             IF (PTEST) THEN
1837:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '1837:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
1838:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1838:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
1839:             END IF1839:             END IF
1840:          ELSE IF (PYT) THEN 
1841:           call py(COORDS,VNEW,ENERGY,.true.) 
1842:             IF (SSTEST) THEN 
1843:                CALL py_secder(COORDS,SSTEST) 
1844:             END IF 
1845:             IF (PTEST) THEN 
1846:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         ' 
1847:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         ' 
1848:             END IF 
1849: 1840: 
1850:          ELSE IF (PTSTSTT) THEN1841:          ELSE IF (PTSTSTT) THEN
1851:             CALL PTSTST(COORDS,VNEW,ENERGY,GTEST,SSTEST)1842:             CALL PTSTST(COORDS,VNEW,ENERGY,GTEST,SSTEST)
1852:             IF (PTEST) THEN1843:             IF (PTEST) THEN
1853:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '1844:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
1854:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1845:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
1855:             END IF1846:             END IF
1856: 1847: 
1857: 1848: 
1858:          ELSE IF (STOCKAAT) THEN1849:          ELSE IF (STOCKAAT) THEN


r31852/rbinertia.f90 2017-01-30 12:30:12.280061097 +0000 r31851/rbinertia.f90 2017-01-30 12:30:16.132113479 +0000
  1:       SUBROUTINE RBINERTIA(Q, ITX, ITY, ITZ)  1:       SUBROUTINE RBINERTIA(Q, ITX, ITY, ITZ)
  2:   2: 
  3:       USE COMMONS, ONLY: NATOMS, NRBSITES, RBSITE  3:       USE COMMONS, ONLY: NATOMS, NRBSITES, RBSITE
  4:       USE KEY, ONLY:     NTSITES, NTIPT, TIPID, STOCKAAT, MULTISITEPYT, PYT  4:       USE KEY, ONLY:     NTSITES, NTIPT, TIPID, STOCKAAT, MULTISITEPYT
  5:   5: 
  6:       IMPLICIT NONE  6:       IMPLICIT NONE
  7:   7: 
  8:       INTEGER          :: J1, J2, J3, J4, J5  8:       INTEGER          :: J1, J2, J3, J4, J5
  9:       DOUBLE PRECISION :: Q(3*NATOMS), ITX, ITY, ITZ, MSRBST(NTSITES), SITEMASS(NRBSITES)  9:       DOUBLE PRECISION :: Q(3*NATOMS), ITX, ITY, ITZ, MSRBST(NTSITES), SITEMASS(NRBSITES)
 10:       DOUBLE PRECISION :: X(3*NTSITES), R(3), P(3), RM(3,3), IT(3,3), CMX, CMY, CMZ, VEC(3,3), MASST 10:       DOUBLE PRECISION :: X(3*NTSITES), R(3), P(3), RM(3,3), IT(3,3), CMX, CMY, CMZ, VEC(3,3), MASST
 11:        11:       
 12:       IF (NTIPT) THEN 12:       IF (NTIPT) THEN
 13:          IF (TIPID == 4) THEN 13:          IF (TIPID == 4) THEN
 14:             SITEMASS(:) = (/16.D0, 1.D0, 1.D0, 0.D0/) 14:             SITEMASS(:) = (/16.D0, 1.D0, 1.D0, 0.D0/)
 15:          ENDIF 15:          ENDIF
 16:       ELSEIF (STOCKAAT.OR.MULTISITEPYT.OR.PYT) THEN 16:       ELSEIF (STOCKAAT.OR.MULTISITEPYT) THEN
 17:          SITEMASS(1) = 1.D0 17:          SITEMASS(1) = 1.D0
 18:       ELSE 18:       ELSE
 19:          SITEMASS(:) = 1.D0  19:          SITEMASS(:) = 1.D0 
 20:       ENDIF 20:       ENDIF
 21:   21:  
 22:       J5 = 0 22:       J5 = 0
 23:  23: 
 24:       DO J1 = 1, NATOMS/2 24:       DO J1 = 1, NATOMS/2
 25:  25: 
 26:          J2   = 3*J1 26:          J2   = 3*J1


r31852/rbperm.f90 2017-01-30 12:30:12.508064197 +0000 r31851/rbperm.f90 2017-01-30 12:30:16.688121156 +0000
 32: !     centre of coordinates of COORDSA will be the same as for COORDSB. 32: !     centre of coordinates of COORDSA will be the same as for COORDSB.
 33: ! 33: !
 34: !     ---------------------------------------------------------------------------------------------- 34: !     ----------------------------------------------------------------------------------------------
 35:   35:  
 36:       SUBROUTINE RBMINPERMDIST(COORDSB,COORDSA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,SITESB,SITESA) 36:       SUBROUTINE RBMINPERMDIST(COORDSB,COORDSA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,SITESB,SITESA)
 37:  37: 
 38: !     returns DISTANCE as the actual distance, rather than the squared distance 38: !     returns DISTANCE as the actual distance, rather than the squared distance
 39:  39: 
 40:       USE COMMONS, ONLY : NATOMS, NRBSITES  40:       USE COMMONS, ONLY : NATOMS, NRBSITES 
 41:       USE KEY,ONLY : NTSITES, NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, GEOMDIFFTOL, BESTPERM, EFIELDT, PAPT, PTSTSTT, & 41:       USE KEY,ONLY : NTSITES, NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, GEOMDIFFTOL, BESTPERM, EFIELDT, PAPT, PTSTSTT, &
 42:      &               RBOPS, NRBGROUP, RBSYMT, DBPTDT, MULTISITEPYT  42:      &               RBOPS, NRBGROUP, RBSYMT, DBPTDT 
 43:  43: 
 44:       IMPLICIT NONE 44:       IMPLICIT NONE
 45:  45: 
 46:       INTEGER, PARAMETER :: MAXIMUMTRIES = 100 46:       INTEGER, PARAMETER :: MAXIMUMTRIES = 100
 47:       INTEGER            :: NPERM, PATOMS, NTRIES, NSIZE, J1, J2, J3, NRB !jwrm2> unused, INFO, JMAX, LOCMAX(1) 47:       INTEGER            :: NPERM, PATOMS, NTRIES, NSIZE, J1, J2, J3, NRB !jwrm2> unused, INFO, JMAX, LOCMAX(1)
 48:       INTEGER            :: INVERT, NORBIT1, NORBIT2, PERM(NATOMS), NCHOOSE2, NDUMMY, LPERM(NATOMS), NCHOOSE1 48:       INTEGER            :: INVERT, NORBIT1, NORBIT2, PERM(NATOMS), NCHOOSE2, NDUMMY, LPERM(NATOMS), NCHOOSE1
 49:       INTEGER            :: NEWPERM(NATOMS), ALLPERM(NATOMS) 49:       INTEGER            :: NEWPERM(NATOMS), ALLPERM(NATOMS)
 50:       DOUBLE PRECISION   :: COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DIST2 !jwrm2> unused, DISTWP, TEMPA(9*NATOMS)   50:       DOUBLE PRECISION   :: COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DIST2 !jwrm2> unused, DISTWP, TEMPA(9*NATOMS)  
 51:       DOUBLE PRECISION   :: DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS) !jwrm2> unused, DUMMYWP(3*NATOMS) 51:       DOUBLE PRECISION   :: DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS) !jwrm2> unused, DUMMYWP(3*NATOMS)
 52:       DOUBLE PRECISION   :: XA(3*NTSITES),  XB(3*NTSITES), XTMP(3*NATOMS) !jwrm2> unused, XBS(3*NTSITES) 52:       DOUBLE PRECISION   :: XA(3*NTSITES),  XB(3*NTSITES), XTMP(3*NATOMS) !jwrm2> unused, XBS(3*NTSITES)
 55:       DOUBLE PRECISION   :: ROTA(3,3), ROTB(3,3) !jwrm2> unused, ROTINVA(3,3), ROTINVB(3,3), RMATCUMUL(3,3),  55:       DOUBLE PRECISION   :: ROTA(3,3), ROTB(3,3) !jwrm2> unused, ROTINVA(3,3), ROTINVB(3,3), RMATCUMUL(3,3), 
 56:       DOUBLE PRECISION   :: RMATBEST(3,3) !jwrm2> unused, RMATWP(3,3), ROTABEST(3,3), ROTINVBBEST(3,3) 56:       DOUBLE PRECISION   :: RMATBEST(3,3) !jwrm2> unused, RMATWP(3,3), ROTABEST(3,3), ROTINVBBEST(3,3)
 57:       DOUBLE PRECISION   :: CMAX, CMAY, CMAZ, CMBX, CMBY, CMBZ 57:       DOUBLE PRECISION   :: CMAX, CMAY, CMAZ, CMBX, CMBY, CMBZ
 58:       DOUBLE PRECISION   :: PDUMMYA(3*NATOMS), PDUMMYB(3*NATOMS), LDISTANCE, XDUMMY, BOXLX, BOXLY, BOXLZ, WORSTRAD 58:       DOUBLE PRECISION   :: PDUMMYA(3*NATOMS), PDUMMYB(3*NATOMS), LDISTANCE, XDUMMY, BOXLX, BOXLY, BOXLZ, WORSTRAD
 59:       DOUBLE PRECISION   :: Q2(4), P(3) !jwrm2> unused, AMAT(4,4), BMAT(4,4), DIAG(4), Q(4), Q1(4) 59:       DOUBLE PRECISION   :: Q2(4), P(3) !jwrm2> unused, AMAT(4,4), BMAT(4,4), DIAG(4), Q(4), Q1(4)
 60:       DOUBLE PRECISION   :: THETAH, DUMMYC(3*NATOMS) !jwrm2> unused, FCT, DUMMYD(3*NATOMS), ST, THETA 60:       DOUBLE PRECISION   :: THETAH, DUMMYC(3*NATOMS) !jwrm2> unused, FCT, DUMMYD(3*NATOMS), ST, THETA
 61:       LOGICAL            :: DEBUG, BULKT 61:       LOGICAL            :: DEBUG, BULKT
 62:       DOUBLE PRECISION   :: BESTA(3*NATOMS), RBDISTANCE, PVEC(3), RTEMP1(3,3), SITEDIST !jwrm2> unused, RTEMP2(3,3) 62:       DOUBLE PRECISION   :: BESTA(3*NATOMS), RBDISTANCE, PVEC(3), RTEMP1(3,3), SITEDIST !jwrm2> unused, RTEMP2(3,3)
 63:       DOUBLE PRECISION   :: QCUMUL(4), QBEST(4), QA(4), QB(4), QBINV(4), QTMP(4) ! jwrm2> unused, QI(4) 63:       DOUBLE PRECISION   :: QCUMUL(4), QBEST(4), QA(4), QB(4), QBINV(4), QTMP(4) ! jwrm2> unused, QI(4)
 64:       DOUBLE PRECISION   :: COORDSBS(3*NATOMS), COORDSAS(3*NATOMS), COORDSCOMA(3*NATOMS/2), COORDSCOMB(3*NATOMS/2) 64:       DOUBLE PRECISION   :: COORDSBS(3*NATOMS), COORDSAS(3*NATOMS), COORDSCOMA(3*NATOMS/2), COORDSCOMB(3*NATOMS/2)
 65:       DOUBLE PRECISION   :: ENERGY, VNEW(3*NATOMS), RMS 
 66:  65: 
 67:       IF (PAPT .OR. PTSTSTT) GOTO 100 66:       IF (PAPT .OR. PTSTSTT) GOTO 100
  67:      
 68:       COORDSCOMA(:) = COORDSA(1:3*NATOMS/2) 68:       COORDSCOMA(:) = COORDSA(1:3*NATOMS/2)
 69:       COORDSCOMB(:) = COORDSB(1:3*NATOMS/2) 69:       COORDSCOMB(:) = COORDSB(1:3*NATOMS/2)
 70:  70: 
 71:       IF (DEBUG) THEN 71:       IF (DEBUG) THEN
 72:          WRITE(*,'(3F20.10)')COORDSA(1:3*NATOMS) 72:          WRITE(*,'(3F20.10)')COORDSA(1:3*NATOMS)
 73:          WRITE(*,*) 73:          WRITE(*,*)
 74:          WRITE(*,'(3F20.10)')COORDSB(1:3*NATOMS) 74:          WRITE(*,'(3F20.10)')COORDSB(1:3*NATOMS)
 75:       ENDIF 75:       ENDIF
 76:  76: 
 77: !     do not call minpermdistrbCOM for multisite PY  77:       NATOMS = NATOMS/2
 78:       IF(MULTISITEPYT) GOTO 100 
 79:         NATOMS = NATOMS/2 
 80:  78: 
 81:         CALL MINPERMDISTRBCOM(COORDSCOMB,COORDSCOMA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT) 79:       CALL  MINPERMDISTRBCOM(COORDSCOMB,COORDSCOMA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT)
 82:  80: 
 83:         NATOMS = 2*NATOMS 81:       NATOMS = 2*NATOMS
 84:  82: 
 85:         IF (SQRT(DISTANCE) <= GEOMDIFFTOL) THEN 83:       IF (SQRT(DISTANCE) <= GEOMDIFFTOL) THEN
 86:          DISTANCE = SQRT(DISTANCE) 84:          DISTANCE = SQRT(DISTANCE)
 87:          IF (DEBUG) PRINT '(A)',' rbpermdist> minpermdistrbcom suggests identical' 85:          IF (DEBUG) PRINT '(A)',' rbpermdist> minpermdistrbcom suggests identical'
 88:          CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0 86:          CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0
 89:          CMBX = 0.0D0; CMBY = 0.0D0; CMBZ = 0.0D0 87:          CMBX = 0.0D0; CMBY = 0.0D0; CMBZ = 0.0D0
 90:          DO J1 = 1, NATOMS/2 88:          DO J1 = 1, NATOMS/2
 91:             CMAX = CMAX + COORDSA(3*(J1-1)+1) 89:             CMAX = CMAX + COORDSA(3*(J1-1)+1)
 92:             CMAY = CMAY + COORDSA(3*(J1-1)+2) 90:             CMAY = CMAY + COORDSA(3*(J1-1)+2)
 93:             CMAZ = CMAZ + COORDSA(3*(J1-1)+3) 91:             CMAZ = CMAZ + COORDSA(3*(J1-1)+3)
 94:             CMBX = CMBX + COORDSB(3*(J1-1)+1) 92:             CMBX = CMBX + COORDSB(3*(J1-1)+1)
 95:             CMBY = CMBY + COORDSB(3*(J1-1)+2) 93:             CMBY = CMBY + COORDSB(3*(J1-1)+2)
104:             COORDSB(3*(J1-1)+1) = COORDSB(3*(J1-1)+1) - CMBX102:             COORDSB(3*(J1-1)+1) = COORDSB(3*(J1-1)+1) - CMBX
105:             COORDSB(3*(J1-1)+2) = COORDSB(3*(J1-1)+2) - CMBY103:             COORDSB(3*(J1-1)+2) = COORDSB(3*(J1-1)+2) - CMBY
106:             COORDSB(3*(J1-1)+3) = COORDSB(3*(J1-1)+3) - CMBZ104:             COORDSB(3*(J1-1)+3) = COORDSB(3*(J1-1)+3) - CMBZ
107:          ENDDO105:          ENDDO
108:          CALL SITEPOS(COORDSB,SITESB)106:          CALL SITEPOS(COORDSB,SITESB)
109:          CALL SITEPOS(COORDSA,SITESA)107:          CALL SITEPOS(COORDSA,SITESA)
110:          DO J1 = 1, NTSITES108:          DO J1 = 1, NTSITES
111:             J2 = 3*J1109:             J2 = 3*J1
112:             SITESA(J2-2:J2) = MATMUL(RMATBEST,SITESA(J2-2:J2))110:             SITESA(J2-2:J2) = MATMUL(RMATBEST,SITESA(J2-2:J2))
113:          ENDDO111:          ENDDO
114:         RETURN112:       RETURN
115:         ENDIF113:       ENDIF
116: 114: 
117: 100   CONTINUE115: 100   CONTINUE
118: 116: 
119:       COORDSBS(1:3*NATOMS) = COORDSB(1:3*NATOMS) ! to trace error, see at the end117:       COORDSBS(1:3*NATOMS) = COORDSB(1:3*NATOMS) ! to trace error, see at the end
120:       COORDSAS(1:3*NATOMS) = COORDSA(1:3*NATOMS) ! to trace error, see at the end118:       COORDSAS(1:3*NATOMS) = COORDSA(1:3*NATOMS) ! to trace error, see at the end
121: !      print *, 'COORDSB'119: !      print *, 'COORDSB'
122: !      print *, COORDSB120: !      print *, COORDSB
123: !      print *, 'COORDSA'121: !      print *, 'COORDSA'
124: !      print *, COORDSA122: !      print *, COORDSA
125: 123: 
126:       NRB   = (NATOMS/2)124:       NRB   = (NATOMS/2)
127:       IF (DBPTDT) THEN125:       IF (DBPTDT) THEN
128:          NSIZE = NTSITES126:          NSIZE = NTSITES
129:       ELSE127:       ELSE
130:          NSIZE = NRB*NRBSITES128:          NSIZE = NRB*NRBSITES
131:       ENDIF129:       ENDIF
132: 130: 
133: !      CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0 
134: !      DO J1 = 1, NATOMS/2 
135: !         CMAX = CMAX + COORDSA(3*(J1-1)+1) 
136: !         CMAY = CMAY + COORDSA(3*(J1-1)+2) 
137: !         CMAZ = CMAZ + COORDSA(3*(J1-1)+3) 
138: !       ENDDO 
139:       CMAX = 2*CMAX/NATOMS; CMAY = 2*CMAY/NATOMS; CMAZ = 2*CMAZ/NATOMS 
140:       CMBX = 0.0D0; CMBY = 0.0D0; CMBZ = 0.0D0131:       CMBX = 0.0D0; CMBY = 0.0D0; CMBZ = 0.0D0
141:       DO J1 = 1, NATOMS/2132:       DO J1 = 1, NATOMS/2
142:          CMBX = CMBX + COORDSB(3*(J1-1)+1)133:          CMBX = CMBX + COORDSB(3*(J1-1)+1)
143:          CMBY = CMBY + COORDSB(3*(J1-1)+2)134:          CMBY = CMBY + COORDSB(3*(J1-1)+2)
144:          CMBZ = CMBZ + COORDSB(3*(J1-1)+3)135:          CMBZ = CMBZ + COORDSB(3*(J1-1)+3)
145:        ENDDO136:        ENDDO
146:       CMBX = 2*CMBX/NATOMS; CMBY = 2*CMBY/NATOMS; CMBZ = 2*CMBZ/NATOMS137:       CMBX = 2*CMBX/NATOMS; CMBY = 2*CMBY/NATOMS; CMBZ = 2*CMBZ/NATOMS
147: ! do not move along the z direction for gravity field (currently implemented for MULTISITEPY) 
148:       IF(EFIELDT.AND.MULTISITEPYT) CMBZ=0.0D0 
149: !138: !
150: !     Bring COORDSB into standard orientation with respect to the site position139: !     Bring COORDSB into standard orientation with respect to the site position
151: !     The standard orientation needs to be done for the sites if we are going to identify140: !     The standard orientation needs to be done for the sites if we are going to identify
152: !     permutation-inversion isomers with respect to the sites metric!141: !     permutation-inversion isomers with respect to the sites metric!
153: !142: !
154: !     DUMMY(1:3*NATOMS)=DUMMYB(1:3*NATOMS)143: !     DUMMY(1:3*NATOMS)=DUMMYB(1:3*NATOMS)
155: !     CALL RBSITESORIENT(DUMMY,DUMMYB,NORBIT1,1,NORBIT2,1,NATOMS,DEBUG,ROTB,ROTINVB)144: !     CALL RBSITESORIENT(DUMMY,DUMMYB,NORBIT1,1,NORBIT2,1,NATOMS,DEBUG,ROTB,ROTINVB)
156: !145: !
157: !     ----------------------------------------------------------------------------------------------146: !     ----------------------------------------------------------------------------------------------
158: !     !147: !     !
159: !     CALL POTENTIAL(COORDSB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)148: !     CALL POTENTIAL(COORDSB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
160: !     PRINT '(2(A,F25.15))','before ENERGYB =',ENERGY,' RMS=',RMS149: !     PRINT '(2(A,F25.15))','before ENERGYB =',ENERGY,' RMS=',RMS
161: !     CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)150: !     CALL POTENTIAL(DUMMYB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
162: !     PRINT '(2(A,F25.15))','before ENERGYA =',ENERGY,' RMS=',RMS 
163: !!     CALL POTENTIAL(DUMMYB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
164: !     PRINT '(2(A,F25.15))',' after ENERGYB =',ENERGY,' RMS=',RMS151: !     PRINT '(2(A,F25.15))',' after ENERGYB =',ENERGY,' RMS=',RMS
165: !     !152: !     !
166: !     ---------------------------------------------------------------------------------------------- 153: !     ---------------------------------------------------------------------------------------------- 
167: !154: !
168: !     CALL SITEPOS(DUMMYB,XB)155: !     CALL SITEPOS(DUMMYB,XB)
169: !     WRITE(975,'(I6)') NRB*NRBSITES156: !     WRITE(975,'(I6)') NRB*NRBSITES
170: !     WRITE(975,'(A,F20.10)') 'B sites 1'157: !     WRITE(975,'(A,F20.10)') 'B sites 1'
171: !     WRITE(975,'(A,3F20.10)') ('LA ',XB(3*(J3-1)+1:3*(J3-1)+3),J3=1,NRB*NRBSITES)158: !     WRITE(975,'(A,3F20.10)') ('LA ',XB(3*(J3-1)+1:3*(J3-1)+3),J3=1,NRB*NRBSITES)
172: 159: 
173: !     Coordinates in DUMMYB do not change but the coordinates in DUMMYA do.160: !     Coordinates in DUMMYB do not change but the coordinates in DUMMYA do.
174: 161: 
175:       DUMMYB(1:3*NATOMS) = COORDSB(1:3*NATOMS)162:       DUMMYB(1:3*NATOMS) = COORDSB(1:3*NATOMS)
176: 163: 
177:       DUMMYC(1:3*NATOMS) = DUMMYB(1:3*NATOMS)164:       DUMMYC(1:3*NATOMS) = DUMMYB(1:3*NATOMS)
178:       CALL RBSITESORIENT(DUMMYC,DUMMY,NORBIT1,1,NORBIT2,1,NATOMS,QB,DEBUG)165:       CALL RBSITESORIENT(DUMMYC,DUMMY,NORBIT1,1,NORBIT2,1,NATOMS,QB,DEBUG)
179:       CALL QROTMAT(QB,ROTB)166:       CALL QROTMAT(QB,ROTB)
180:       DUMMYB(1:3*NATOMS)= DUMMY(1:3*NATOMS)167:       DUMMYB(1:3*NATOMS)= DUMMY(1:3*NATOMS)
181:       CALL SITEPOS(DUMMYB,XB)168:       CALL SITEPOS(DUMMYB,XB)
182: !169: 
183:       INVERT = 1170:       INVERT = 1
184: 171: 
185:       DBEST    = 1.0D100172:       DBEST    = 1.0D100
186: 60    NCHOOSE1 = 0173: 60    NCHOOSE1 = 0
187: 65    NCHOOSE1 = NCHOOSE1+1174: 65    NCHOOSE1 = NCHOOSE1+1
188: 40    NCHOOSE2 = 0175: 40    NCHOOSE2 = 0
189: 30    NCHOOSE2 = NCHOOSE2+1176: 30    NCHOOSE2 = NCHOOSE2+1
190: 177: 
191:       DO J1 = 1, NATOMS178:       DO J1 = 1, NATOMS
192:          ALLPERM(J1) = J1179:          ALLPERM(J1) = J1
229: !         DISTANCE = DISTANCE + (DUMMYA(J1) - DUMMYB(J1))**2216: !         DISTANCE = DISTANCE + (DUMMYA(J1) - DUMMYB(J1))**2
230: !      ENDDO217: !      ENDDO
231: !     !218: !     !
232: !     ----------------------------------------------------------------------------------------------219: !     ----------------------------------------------------------------------------------------------
233:       220:       
234: !      IF (INVERT == -1) ROTA(:,:) = MATMUL(ROTA(:,:),RMATI(:,:))  221: !      IF (INVERT == -1) ROTA(:,:) = MATMUL(ROTA(:,:),RMATI(:,:))  
235: !     RMATI was returned on inversion followed by rotation to bring about best superposition on DUMMYB in one go222: !     RMATI was returned on inversion followed by rotation to bring about best superposition on DUMMYB in one go
236: 223: 
237:       DUMMYA(1:3*NATOMS) = COORDSA(1:3*NATOMS)224:       DUMMYA(1:3*NATOMS) = COORDSA(1:3*NATOMS)
238:       DUMMYC(1:3*NATOMS) = DUMMYA(1:3*NATOMS)225:       DUMMYC(1:3*NATOMS) = DUMMYA(1:3*NATOMS)
239:  
240:       CALL RBSITESORIENT(DUMMYC,DUMMY,NORBIT1,NCHOOSE1,NORBIT2,NCHOOSE2,NATOMS,QA,DEBUG)226:       CALL RBSITESORIENT(DUMMYC,DUMMY,NORBIT1,NCHOOSE1,NORBIT2,NCHOOSE2,NATOMS,QA,DEBUG)
241:       CALL QROTMAT(QA,ROTA)227:       CALL QROTMAT(QA,ROTA)
242:       DUMMYA(1:3*NATOMS) = DUMMY(1:3*NATOMS)228:       DUMMYA(1:3*NATOMS) = DUMMY(1:3*NATOMS)
243:       CALL SITEPOS(DUMMYA,XA)229:       CALL SITEPOS(DUMMYA,XA)
244: 230: 
245: !      IF (INVERT == -1) THEN231: !      IF (INVERT == -1) THEN
246: !         QTMP(1:4) = QI(1:4)232: !         QTMP(1:4) = QI(1:4)
247: !         CALL QROTQ(QA,QTMP)233: !         CALL QROTQ(QA,QTMP)
248: !         QA(1:4) = QTMP(1:4)234: !         QA(1:4) = QTMP(1:4)
249: !      ENDIF  235: !      ENDIF  
452:                DUMMYA(1:6*NRB)=BESTA(1:6*NRB)438:                DUMMYA(1:6*NRB)=BESTA(1:6*NRB)
453:             ENDDO439:             ENDDO
454: !440: !
455: !     This call aligns the overall orientation with respect to the sites metric. We should already441: !     This call aligns the overall orientation with respect to the sites metric. We should already
456: !     have identified permutation-inversion isomers using the alignment based on centres of mass442: !     have identified permutation-inversion isomers using the alignment based on centres of mass
457: !     above, if applicable. We now repeat but with an orientational alignment with respect to sites.443: !     above, if applicable. We now repeat but with an orientational alignment with respect to sites.
458: !444: !
459:             CALL RBMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,Q2,DEBUG)445:             CALL RBMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,Q2,DEBUG)
460:  446:  
461:             CALL QROTQ(Q2,QCUMUL)447:             CALL QROTQ(Q2,QCUMUL)
 448: 
462:             DO J1=1,NRB449:             DO J1=1,NRB
463:                DO J2=1,NRBGROUP450:                DO J2=1,NRBGROUP
464: !451: !
465: !     MUST NOT USE RMAT HERE - because RMAT is accumulated below.452: !     MUST NOT USE RMAT HERE - because RMAT is accumulated below.
466: ! 453: ! 
467:                   P(:)      = DUMMYA(3*(NRB+J1-1)+1:3*(NRB+J1-1)+3)454:                   P(:)      = DUMMYA(3*(NRB+J1-1)+1:3*(NRB+J1-1)+3)
468:                   IF (J2 /= 1) THEN455:                   IF (J2 /= 1) THEN
469:                      CALL ROTMAT(P,RTEMP1)456:                      CALL ROTMAT(P,RTEMP1)
470:                      PVEC(1:3) = MATMUL(RTEMP1,RBOPS(1:3,J2))457:                      PVEC(1:3) = MATMUL(RTEMP1,RBOPS(1:3,J2))
471:                      PVEC(1:3) = PVEC(1:3)/SQRT(DOT_PRODUCT(PVEC(1:3),PVEC(1:3)))458:                      PVEC(1:3) = PVEC(1:3)/SQRT(DOT_PRODUCT(PVEC(1:3),PVEC(1:3)))
479: !466: !
480: !                 CALL POTENTIAL(DUMMYC,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)467: !                 CALL POTENTIAL(DUMMYC,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
481: !                 PRINT '(A,I6,2(A,F25.15))',' rbminpermdist> after RBROT following internal symmetry 468: !                 PRINT '(A,I6,2(A,F25.15))',' rbminpermdist> after RBROT following internal symmetry 
482: !                                              operation ',J2,' energy for A=', ENERGY,' RMS=',RMS469: !                                              operation ',J2,' energy for A=', ENERGY,' RMS=',RMS
483: !470: !
484: !     Calculate new distance based on all sites metric rather than just centre-of-mass rigid body coordinates.471: !     Calculate new distance based on all sites metric rather than just centre-of-mass rigid body coordinates.
485: ! 472: ! 
486:                   CALL SITEPOS(DUMMYC,SITESA)473:                   CALL SITEPOS(DUMMYC,SITESA)
487:                   SITEDIST=0.D0474:                   SITEDIST=0.D0
488:                   DO J3=1,3*NRB*NRBSITES475:                   DO J3=1,3*NRB*NRBSITES
489:                        SITEDIST=SITEDIST+(SITESA(J3)-SITESB(J3))**2476:                      SITEDIST=SITEDIST+(SITESA(J3)-SITESB(J3))**2
490:                     ENDDO477:                   ENDDO
491:   !                 PRINT '(A,2F20.10)',' rbminpermdist> sites distance^2,RBDISTANCE=',SITEDIST,RBDISTANCE478: !                 PRINT '(A,2F20.10)',' rbminpermdist> sites distance^2,RBDISTANCE=',SITEDIST,RBDISTANCE
492:                     IF (SITEDIST.LT.RBDISTANCE) THEN479:                   IF (SITEDIST.LT.RBDISTANCE) THEN
493:                        RBDISTANCE=SITEDIST480:                      RBDISTANCE=SITEDIST
494:                        BESTA(1:6*NRB)=DUMMYC(1:6*NRB)481:                      BESTA(1:6*NRB)=DUMMYC(1:6*NRB)
495:                     ENDIF482:                   ENDIF
496:                  ENDDO483:                ENDDO
497:                  DUMMYA(1:6*NRB)=BESTA(1:6*NRB)484:                DUMMYA(1:6*NRB)=BESTA(1:6*NRB)
498:               ENDDO485:             ENDDO
499:            ELSE486:          ELSE
500:   !487: !
501:   !     This call aligns the overall orientation with respect to the sites metric.488: !     This call aligns the overall orientation with respect to the sites metric.
502:   !     Internal symmetry operations of the rigid bodies are not considered.489: !     Internal symmetry operations of the rigid bodies are not considered.
503:   !490: !
504:   !sf344 debug491:             CALL RBMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,Q2,DEBUG)
505:   !     CALL POTENTIAL(DUMMYB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)492: !     accumulate rotation
506:   !     PRINT '(2(A,F25.15))','DUMMYB before RBMINDIST =',ENERGY,' RMS=',RMS493:             CALL QROTQ(Q2,QCUMUL)
507:   !     CALL POTENTIAL(DUMMYA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)494: 
508:   !     PRINT '(2(A,F25.15))','DUMMYA before RBMINDIST 1 =',ENERGY,' RMS=',RMS495:          ENDIF
509: 496: 
510:               CALL RBMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,Q2,DEBUG)497:          DISTANCE  = DISTANCE*DISTANCE
511: 498: 
512:   !     accumulate rotation499:          IF (NTRIES .LT. MAXIMUMTRIES) THEN
513:               CALL QROTQ(Q2,QCUMUL)500:             GOTO 10
514: 501: !         ELSE ! prevent infinite loop
515:            ENDIF502: !            IF (INVERT == -1) THEN
516: 503: !               IF(DEBUG) PRINT '(A)',' rbminpermdist> WARNING - number of tries exceeded, giving up'
517:            DISTANCE  = DISTANCE*DISTANCE504: !            ELSE
518: 505: !               PRINT '(A)',' rbminpermdist> WARNING - number of tries exceeded, giving up'
519:            IF (NTRIES .LT. MAXIMUMTRIES) THEN506: !            ENDIF
520:               GOTO 10507:          ENDIF
521:   !         ELSE ! prevent infinite loop508: 
522:   !            IF (INVERT == -1) THEN509:       ENDIF
523:   !               IF(DEBUG) PRINT '(A)',' rbminpermdist> WARNING - number of tries exceeded, giving up'510: 
524:   !            ELSE511:       IF (DISTANCE .LT. DBEST) THEN
525:   !               PRINT '(A)',' rbminpermdist> WARNING - number of tries exceeded, giving up'512: 
526:   !            ENDIF513:          DBEST                = DISTANCE
527:            ENDIF514:          XBEST(1:3*NATOMS)    = DUMMYA(1:3*NATOMS)
528: 515:          BESTPERM(1:NATOMS)   = ALLPERM(1:NATOMS)
529:         ENDIF516:          QTMP(1:4)  = QA(1:4)
530: 517:          CALL QROTQ(QCUMUL,QTMP)
531:         IF (DISTANCE .LT. DBEST) THEN518:          QBEST(1:4) = QTMP(1:4)
532: 519:          QBINV(1:4) = (/QB(1), -QB(2:4)/)
533:            DBEST                = DISTANCE520: 
534:            XBEST(1:3*NATOMS)    = DUMMYA(1:3*NATOMS)521:       ENDIF
535:            BESTPERM(1:NATOMS)   = ALLPERM(1:NATOMS)522: !
536:            QTMP(1:4)  = QA(1:4)523: ! If GEOMDIFFTOLL is set too sloppy we could miss the best solution by exiting via the line
537:            CALL QROTQ(QCUMUL,QTMP)524: ! below without trying other orbits. Turn off this escape?!
538:            QBEST(1:4) = QTMP(1:4)525: ! The turn off seems to be a bug for non-RBPERM runs. LJ38 tests fail with all compilers, producing
539:            QBINV(1:4) = (/QB(1), -QB(2:4)/)526: ! a "Distance is zero: this should not happen" error!!! DJW
540: 527: ! Put the escape back in for now.
541:         ENDIF528: 
542:   !529:       IF (SQRT(DBEST).LT.GEOMDIFFTOL) GOTO 50
543:   ! If GEOMDIFFTOLL is set too sloppy we could miss the best solution by exiting via the line530:       IF (NCHOOSE2.LT.NORBIT2) GOTO 30
544:   ! below without trying other orbits. Turn off this escape?!531:       IF (NCHOOSE1.LT.NORBIT1) GOTO 65
545:   ! The turn off seems to be a bug for non-RBPERM runs. LJ38 tests fail with all compilers, producing532:       IF (EFIELDT) GOTO 50
546:   ! a "Distance is zero: this should not happen" error!!! DJW533: !
547:   ! Put the escape back in for now.534: 50    DISTANCE = DBEST       ! squared distance
548: 535: !
549:         IF (SQRT(DBEST).LT.GEOMDIFFTOL) GOTO 50536: !     XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.
550:         IF (NCHOOSE2.LT.NORBIT2) GOTO 30537: !     Rotate XBEST by ROTINVBBEST to put in best correspondence with COORDSB, undoing the reorientation to  
551:         IF (NCHOOSE1.LT.NORBIT1) GOTO 65538: !     DUMMYB from RBORIENT. 
552:         IF (EFIELDT) GOTO 50539: !     We should get the same result for ROTINVBBEST * RMATBEST * (COORDSA-CMA), 
553:   !540: !     where RMATBEST = +/- RMATCUMUL * ROTA for the best alignment 
554:   50    DISTANCE = DBEST       ! squared distance541: !     (aside from a possible permutation of the atom ordering)
555:   !542:       CALL RBROT(XBEST, XTMP, QBINV, NATOMS)
556:   !     XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.543:       XBEST(1:3*NATOMS) = XTMP(1:3*NATOMS)
557:   !     Rotate XBEST by ROTINVBBEST to put in best correspondence with COORDSB, undoing the reorientation to  544: 
558:   !     DUMMYB from RBORIENT. 545:       DO J1 = 1, (NATOMS/2)
559:   !     We should get the same result for ROTINVBBEST * RMATBEST * (COORDSA-CMA), 546:          XBEST(3*(J1-1)+1) = XBEST(3*(J1-1)+1) + CMBX
560:   !     where RMATBEST = +/- RMATCUMUL * ROTA for the best alignment 547:          XBEST(3*(J1-1)+2) = XBEST(3*(J1-1)+2) + CMBY
561:   !     (aside from a possible permutation of the atom ordering)548:          XBEST(3*(J1-1)+3) = XBEST(3*(J1-1)+3) + CMBZ
562: 549:       ENDDO
563:         CALL RBROT(XBEST, XTMP, QBINV, NATOMS)550: 
564:         XBEST(1:3*NATOMS) = XTMP(1:3*NATOMS)551: !     Site-site distance
565: 552: 
566:         DO J1 = 1, (NATOMS/2)553:       CALL SITEPOS(XBEST,XA)
567:            XBEST(3*(J1-1)+1) = XBEST(3*(J1-1)+1) + CMBX554: 
568:            XBEST(3*(J1-1)+2) = XBEST(3*(J1-1)+2) + CMBY555:       CALL SITEPOS(COORDSB,XB)
569:            XBEST(3*(J1-1)+3) = XBEST(3*(J1-1)+3) + CMBZ556: 
570:         ENDDO557:       XDUMMY = 0.D0
571: 558:       DO J1 = 1, 3*NSIZE
572: !       CALL POTENTIAL(XBEST,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)559:          XDUMMY = XDUMMY + (XA(J1) - XB(J1))**2
573:    560:       ENDDO
574:   !     Site-site distance561: 
575: 562:       IF (ABS(SQRT(XDUMMY)-SQRT(DISTANCE)).GT.GEOMDIFFTOL) THEN
576:         CALL SITEPOS(XBEST,XA)563:           PRINT '(2(A,F20.10))','rbminpermdist> ERROR *** distance between transformed XBEST and COORDSB=',  &
577: 564:      &    SQRT(XDUMMY),  ' should be ', SQRT(DISTANCE)
578:         CALL SITEPOS(COORDSB,XB)565:           PRINT *, 'COORDSB'
579: 566:           DO J1 = 1, NATOMS
580:         XDUMMY = 0.D0567:              J2 = 3*J1
581:         DO J1 = 1, 3*NSIZE568:              PRINT *, COORDSBS(J2-2), COORDSBS(J2-1), COORDSBS(J2)
582:            XDUMMY = XDUMMY + (XA(J1) - XB(J1))**2569:           ENDDO 
583:         ENDDO570:           PRINT *, 'COORDSA'
584: 571:           DO J1 = 1, NATOMS
585:         IF (ABS(SQRT(XDUMMY)-SQRT(DISTANCE)).GT.GEOMDIFFTOL) THEN572:              J2 = 3*J1
586:             PRINT '(2(A,F20.10))','rbminpermdist> ERROR *** distance between transformed XBEST and COORDSB=',  &573:              PRINT *, COORDSAS(J2-2), COORDSAS(J2-1), COORDSAS(J2)
587:        &    SQRT(XDUMMY),  ' should be ', SQRT(DISTANCE)574:           ENDDO 
588:             PRINT *, 'COORDSB'575: !          PRINT *, 'XBEST'
589:             DO J1 = 1, NATOMS576: !          DO J1 = 1, NATOMS
590:                J2 = 3*J1577: !             J2 = 3*J1
591:                PRINT *, COORDSBS(J2-2), COORDSBS(J2-1), COORDSBS(J2)578: !             PRINT *, XBEST(J2-2), XBEST(J2-1), XBEST(J2)
592:             ENDDO 579: !          ENDDO 
593:             PRINT *, 'COORDSA'580: 
594:             DO J1 = 1, NATOMS581: !          PRINT *, 'RMSD=', SQRT(DISTANCE/DFLOAT(NATOMS*NRBSITES/2))
595:                J2 = 3*J1582:          STOP ! commented to allow newtip to continue as a temporary workaround DJW 2/12/12
596:                PRINT *, COORDSAS(J2-2), COORDSAS(J2-1), COORDSAS(J2)583: 
597:             ENDDO 584:       ENDIF
598:   !          PRINT *, 'XBEST'585: 
599:   !          DO J1 = 1, NATOMS586:       CALL QROTQ(QBINV,QBEST)
600:   !             J2 = 3*J1587: 
601:   !             PRINT *, XBEST(J2-2), XBEST(J2-1), XBEST(J2)588:       CALL QROTMAT(QBEST,RMATBEST)
602:   !          ENDDO 589: 
603: 590:       COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!
604:   !          PRINT *, 'RMSD=', SQRT(DISTANCE/DFLOAT(NATOMS*NRBSITES/2))591: 
605:            STOP ! commented to allow newtip to continue as a temporary workaround DJW 2/12/12592: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!debug
606: 593: !      CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
607:         ENDIF594: !      PRINT '(2(A,F25.15))',' finally A energy=',ENERGY,' RMS=',RMS
608:  
609:         CALL QROTQ(QBINV,QBEST) 
610:  
611:         CALL QROTMAT(QBEST,RMATBEST) 
612:  
613:         COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input! 
614:  
615:   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!debug 
616: !        CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
617: !        PRINT '(2(A,F25.15))',' finally A energy=',ENERGY,' RMS=',RMS 
618: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!debug595: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!debug
619: 596: 
620:       DISTANCE=SQRT(DISTANCE)597:       DISTANCE=SQRT(DISTANCE)
621: 598: 
622:       IF (DISTANCE <= GEOMDIFFTOL .AND. (PAPT .OR. PTSTSTT)) THEN599:       IF (DISTANCE <= GEOMDIFFTOL .AND. (PAPT .OR. PTSTSTT)) THEN
623:          CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0600:          CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0
624:          CMBX = 0.0D0; CMBY = 0.0D0; CMBZ = 0.0D0601:          CMBX = 0.0D0; CMBY = 0.0D0; CMBZ = 0.0D0
625:          DO J1 = 1, NATOMS/2602:          DO J1 = 1, NATOMS/2
626:             CMAX = CMAX + COORDSAS(3*(J1-1)+1)603:             CMAX = CMAX + COORDSAS(3*(J1-1)+1)
627:             CMAY = CMAY + COORDSAS(3*(J1-1)+2)604:             CMAY = CMAY + COORDSAS(3*(J1-1)+2)
695:             ENDDO672:             ENDDO
696:          ENDIF673:          ENDIF
697:       ENDDO674:       ENDDO
698:      675:      
699:       END SUBROUTINE SITEPOS676:       END SUBROUTINE SITEPOS
700: 677: 
701: !     ----------------------------------------------------------------------------------------------678: !     ----------------------------------------------------------------------------------------------
702: 679: 
703:       SUBROUTINE RBSITESORIENT(X, T1, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2, NATOMS, Q2, DEBUG)680:       SUBROUTINE RBSITESORIENT(X, T1, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2, NATOMS, Q2, DEBUG)
704: 681: 
705:       USE KEY, ONLY : NTSITES, EFIELDT, MULTISITEPYT 682:       USE KEY, ONLY : NTSITES, EFIELDT 
706: !                        683: !                        
707: !     This subroutine puts permutational isomers of a rigid-body system into a standard orientation684: !     This subroutine puts permutational isomers of a rigid-body system into a standard orientation
708: !     with respect to the sites685: !     with respect to the sites
709: !686: !
710:       IMPLICIT NONE687:       IMPLICIT NONE
711:       INTEGER          :: NATOMS, I, J, J1, J2, JMAX1, JMAX2, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2, NSIZE !jwrm2> unused, OFFSET688:       INTEGER          :: NATOMS, I, J, J1, J2, JMAX1, JMAX2, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2, NSIZE !jwrm2> unused, OFFSET
712:       DOUBLE PRECISION :: X(3*NATOMS), T1(3*NATOMS)689:       DOUBLE PRECISION :: X(3*NATOMS), T1(3*NATOMS)
713: !      DOUBLE PRECISION :: XS(3*NATOMS*NRBSITES/2), T1S(3*NATOMS*NRBSITES/2), T2S(3*NATOMS*NRBSITES/2), DIST(NATOMS*NRBSITES/2)690: !      DOUBLE PRECISION :: XS(3*NATOMS*NRBSITES/2), T1S(3*NATOMS*NRBSITES/2), T2S(3*NATOMS*NRBSITES/2), DIST(NATOMS*NRBSITES/2)
714:       DOUBLE PRECISION :: XS(3*NTSITES), T1S(3*NTSITES), T2S(3*NTSITES), DIST(NTSITES)691:       DOUBLE PRECISION :: XS(3*NTSITES), T1S(3*NTSITES), T2S(3*NTSITES), DIST(NTSITES)
715:       DOUBLE PRECISION :: AX(3), Q2(4), ROTM(3,3) !jwrm2> unused, P(3), ROTMINV(3,3)692:       DOUBLE PRECISION :: AX(3), Q2(4), ROTM(3,3) !jwrm2> unused, P(3), ROTMINV(3,3)
716:       DOUBLE PRECISION :: THETA, THETAH, FCT !jwrm2> unused, COST, COSTH, SINT, SINTH, ST693:       DOUBLE PRECISION :: THETA, THETAH, FCT !jwrm2> unused, COST, COSTH, SINT, SINTH, ST
717:       DOUBLE PRECISION :: CMX, CMY, CMZ, DMAX, PROJ, DMAX2, CUTOFF1, DTEMP !jwrm2> unused, DUMMY694:       DOUBLE PRECISION :: CMX, CMY, CMZ, DMAX, PROJ, DMAX2, CUTOFF1, DTEMP !jwrm2> unused, DUMMY
718:       LOGICAL          :: DEBUG695:       LOGICAL          :: DEBUG
719:       DOUBLE PRECISION :: ENERGY, VNEW(3*NATOMS), RMS 
720: 696: 
721:       NSIZE = NTSITES697:       NSIZE = NTSITES
722:       !jwrm2> unused: OFFSET  = 3*NATOMS/2698:       !jwrm2> unused: OFFSET  = 3*NATOMS/2
723:       CUTOFF1 = 1.D-03699:       CUTOFF1 = 1.D-03
724: 700: 
725:       CALL SITEPOS(X, XS)701:       CALL SITEPOS(X, XS)
726: !702: !
727: !     Move centre of mass to the origin.703: !     Move centre of mass to the origin.
728: !704: !
729:       CMX = 0.D0; CMY = 0.D0; CMZ = 0.D0705:       CMX = 0.D0; CMY = 0.D0; CMZ = 0.D0
732:          CMX = CMX + XS(J-2)708:          CMX = CMX + XS(J-2)
733:          CMY = CMY + XS(J-1)709:          CMY = CMY + XS(J-1)
734:          CMZ = CMZ + XS(J)710:          CMZ = CMZ + XS(J)
735:       ENDDO711:       ENDDO
736: 712: 
737:       CMX = CMX/NSIZE; CMY = CMY/NSIZE; CMZ = CMZ/NSIZE713:       CMX = CMX/NSIZE; CMY = CMY/NSIZE; CMZ = CMZ/NSIZE
738:       DO I = 1, NSIZE714:       DO I = 1, NSIZE
739:          J = 3*I715:          J = 3*I
740:          XS(J-2) = XS(J-2) - CMX716:          XS(J-2) = XS(J-2) - CMX
741:          XS(J-1) = XS(J-1) - CMY717:          XS(J-1) = XS(J-1) - CMY
742: ! sf344> modification 
743: !        IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)  
744:          XS(J)   = XS(J)   - CMZ718:          XS(J)   = XS(J)   - CMZ
745:       ENDDO719:       ENDDO
746: 720: 
747:       DMAX    = -1.D0721:       DMAX    = -1.D0
748:       NORBIT1 = 1722:       NORBIT1 = 1
749: 723: 
750:       IF (EFIELDT) THEN724:       IF (EFIELDT) THEN
751:          T1S(:) = XS(:)725:          T1S(:) = XS(:)
752:          
753:          T1(:) = X(:) 
754:          GOTO 100726:          GOTO 100
755:       ENDIF727:       ENDIF
756: !728: !
757: !     Find the rbsite which is at the largest distance from the centre of mass729: !     Find the rbsite which is at the largest distance from the centre of mass
758: !730: !
759:       DO J1 = 1, NSIZE731:       DO J1 = 1, NSIZE
760: 732: 
761:          J = 3*J1733:          J = 3*J1
762:          DIST(J1) = SQRT(XS(J-2)**2 + XS(J-1)**2 + XS(J)**2)734:          DIST(J1) = SQRT(XS(J-2)**2 + XS(J-1)**2 + XS(J)**2)
763: 735: 
989: 961: 
990:       SINV(:,:) = A(:,:)/DET962:       SINV(:,:) = A(:,:)/DET
991: 963: 
992:       END SUBROUTINE INVMTRX964:       END SUBROUTINE INVMTRX
993: 965: 
994: !     ----------------------------------------------------------------------------------------------966: !     ----------------------------------------------------------------------------------------------
995: 967: 
996:       SUBROUTINE RBROT(T, X, Q2, NATOMS)968:       SUBROUTINE RBROT(T, X, Q2, NATOMS)
997: !     takes the set of rigid-body coordinates T and returns X after rotation via the quaternion Q2 969: !     takes the set of rigid-body coordinates T and returns X after rotation via the quaternion Q2 
998: !     about the origin970: !     about the origin
999:       USE KEY, ONLY : EFIELDT, MULTISITEPYT971: 
1000:       IMPLICIT NONE972:       IMPLICIT NONE
1001: 973: 
1002:       INTEGER          :: I, J, NATOMS974:       INTEGER          :: I, J, NATOMS
1003:       DOUBLE PRECISION :: T(3*NATOMS), X(3*NATOMS), T1(1:3), Q2(4), P(3), RM(3,3) 975:       DOUBLE PRECISION :: T(3*NATOMS), X(3*NATOMS), T1(1:3), Q2(4), P(3), RM(3,3) 
1004:       DOUBLE PRECISION :: CMX, CMY, CMZ976:       DOUBLE PRECISION :: CMX, CMY, CMZ
1005: !977: !
1006: !     Move centre of mass to the origin.978: !     Move centre of mass to the origin.
1007: !979: !
1008:       CMX = 0.D0; CMY = 0.D0; CMZ = 0.D0980:       CMX = 0.D0; CMY = 0.D0; CMZ = 0.D0
1009:       DO I = 1, NATOMS/2981:       DO I = 1, NATOMS/2
1010:          J = 3*I982:          J = 3*I
1011:          CMX = CMX + T(J-2)983:          CMX = CMX + T(J-2)
1012:          CMY = CMY + T(J-1)984:          CMY = CMY + T(J-1)
1013:          CMZ = CMZ + T(J)985:          CMZ = CMZ + T(J)
1014:       ENDDO986:       ENDDO
1015:       CMX = 2*CMX/NATOMS; CMY = 2*CMY/NATOMS; CMZ = 2*CMZ/NATOMS987:       CMX = 2*CMX/NATOMS; CMY = 2*CMY/NATOMS; CMZ = 2*CMZ/NATOMS
1016:       DO I = 1, NATOMS/2988:       DO I = 1, NATOMS/2
1017:          J      = 3*I989:          J      = 3*I
1018:          X(J-2) = T(J-2) - CMX990:          X(J-2) = T(J-2) - CMX
1019:          X(J-1) = T(J-1) - CMY991:          X(J-1) = T(J-1) - CMY
1020:          IF(EFIELDT.AND.MULTISITEPYT) CMZ = 0.0D0  
1021:          X(J)   = T(J)   - CMZ992:          X(J)   = T(J)   - CMZ
1022:       ENDDO993:       ENDDO
1023: 994: 
1024: !     extract the rotation matrix corresponding to rotation via Q995: !     extract the rotation matrix corresponding to rotation via Q
1025: 996: 
1026:       CALL QROTMAT(Q2,RM)997:       CALL QROTMAT(Q2,RM)
1027: 998: 
1028:       DO I = 1, NATOMS/2999:       DO I = 1, NATOMS/2
1029: 1000: 
1030:          J        = 3*I1001:          J        = 3*I
1120:       RM(3,1) = 2.D0*(Q(2)*Q(4) - Q(1)*Q(3))1091:       RM(3,1) = 2.D0*(Q(2)*Q(4) - Q(1)*Q(3))
1121:       RM(3,2) = 2.D0*(Q(3)*Q(4) + Q(1)*Q(2))1092:       RM(3,2) = 2.D0*(Q(3)*Q(4) + Q(1)*Q(2))
1122:       RM(3,3) = Q(1)**2 + Q(4)**2 - Q(2)**2 - Q(3)**21093:       RM(3,3) = Q(1)**2 + Q(4)**2 - Q(2)**2 - Q(3)**2
1123: 1094: 
1124:       END SUBROUTINE1095:       END SUBROUTINE
1125: 1096: 
1126: !     ----------------------------------------------------------------------------------------------1097: !     ----------------------------------------------------------------------------------------------
1127: 1098: 
1128:       SUBROUTINE RBORIENT(X, T1, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2, NATOMS, Q2, DEBUG)1099:       SUBROUTINE RBORIENT(X, T1, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2, NATOMS, Q2, DEBUG)
1129:       1100:       
1130:       USE KEY, ONLY : EFIELDT,MULTISITEPYT 1101:       USE KEY, ONLY : EFIELDT 
1131: !                        1102: !                        
1132: !     This subroutine puts permutational isomers of a rigid-body system into a standard orientation.1103: !     This subroutine puts permutational isomers of a rigid-body system into a standard orientation.
1133: !1104: !
1134:       IMPLICIT NONE1105:       IMPLICIT NONE
1135:       INTEGER          :: NATOMS, I, J, J1, J2, JMAX1, JMAX2, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2 !jwrm2> unused, OFFSET1106:       INTEGER          :: NATOMS, I, J, J1, J2, JMAX1, JMAX2, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2 !jwrm2> unused, OFFSET
1136:       DOUBLE PRECISION :: X(3*NATOMS), DIST(NATOMS), T1(3*NATOMS), T2(3*NATOMS)1107:       DOUBLE PRECISION :: X(3*NATOMS), DIST(NATOMS), T1(3*NATOMS), T2(3*NATOMS)
1137:       DOUBLE PRECISION :: AX(3), Q2(4), ROTM(3,3) !jwrm2> unused, P(3), Q(4), Q1(4), RM(3,3), ROTMINV(3,3)1108:       DOUBLE PRECISION :: AX(3), Q2(4), ROTM(3,3) !jwrm2> unused, P(3), Q(4), Q1(4), RM(3,3), ROTMINV(3,3)
1138:       DOUBLE PRECISION :: THETA, THETAH, FCT !jwrm2> unused, COST, COSTH, SINT, SINTH, ST1109:       DOUBLE PRECISION :: THETA, THETAH, FCT !jwrm2> unused, COST, COSTH, SINT, SINTH, ST
1139:       DOUBLE PRECISION :: CMX, CMY, CMZ, DMAX, PROJ, DMAX2, CUTOFF1, DTEMP !jwrm2> unused, DUMMY1110:       DOUBLE PRECISION :: CMX, CMY, CMZ, DMAX, PROJ, DMAX2, CUTOFF1, DTEMP !jwrm2> unused, DUMMY
1140:       LOGICAL          :: DEBUG1111:       LOGICAL          :: DEBUG
1150:          CMX = CMX + X(J-2)1121:          CMX = CMX + X(J-2)
1151:          CMY = CMY + X(J-1)1122:          CMY = CMY + X(J-1)
1152:          CMZ = CMZ + X(J)1123:          CMZ = CMZ + X(J)
1153:       ENDDO1124:       ENDDO
1154: 1125: 
1155:       CMX = 2*CMX/NATOMS; CMY = 2*CMY/NATOMS; CMZ = 2*CMZ/NATOMS1126:       CMX = 2*CMX/NATOMS; CMY = 2*CMY/NATOMS; CMZ = 2*CMZ/NATOMS
1156:       DO I = 1, NATOMS/21127:       DO I = 1, NATOMS/2
1157:          J = 3*I1128:          J = 3*I
1158:          X(J-2) = X(J-2) - CMX1129:          X(J-2) = X(J-2) - CMX
1159:          X(J-1) = X(J-1) - CMY1130:          X(J-1) = X(J-1) - CMY
1160: ! do not move along the z direction for gravity field (currently implemented in MULTISITEPY)1131:          X(J)   = X(J)   - CMZ
1161:         IF(EFIELDT.AND.MULTISITEPYT) CMZ=0.0D0 
1162:         X(J)   = X(J)   - CMZ 
1163:       ENDDO1132:       ENDDO
1164: 1133: 
1165:       DMAX    = -1.D01134:       DMAX    = -1.D0
1166:       NORBIT1 = 11135:       NORBIT1 = 1
1167: 1136: 
1168:       IF (EFIELDT) THEN1137:       IF (EFIELDT) THEN
1169: 1138: 
1170:          T1(:) = X(:)1139:          T1(:) = X(:)
1171:          GOTO 1001140:          GOTO 100
1172: 1141: 
1382:       T1(1:3*NATOMS) = X(1:3*NATOMS)1351:       T1(1:3*NATOMS) = X(1:3*NATOMS)
1383: 1352: 
1384:       END SUBROUTINE RBROTXZ1353:       END SUBROUTINE RBROTXZ
1385: 1354: 
1386: !     ----------------------------------------------------------------------------------------------1355: !     ----------------------------------------------------------------------------------------------
1387: 1356: 
1388:       SUBROUTINE RBFRAME(C,XS,Q,RMATBEST)1357:       SUBROUTINE RBFRAME(C,XS,Q,RMATBEST)
1389: 1358: 
1390:       USE COMMONS, ONLY: NATOMS, NRBSITES1359:       USE COMMONS, ONLY: NATOMS, NRBSITES
1391:       USE KEY, ONLY    : DBPT, DBPTDT, DMBLPYT, NTIPT, MSSTOCKT, PAHAT, PATCHYDT, PAPT, PTSTSTT, NTSITES, NCARBON, &1360:       USE KEY, ONLY    : DBPT, DBPTDT, DMBLPYT, NTIPT, MSSTOCKT, PAHAT, PATCHYDT, PAPT, PTSTSTT, NTSITES, NCARBON, &
1392: &                        NPYSITE, MULTISITEPYT, PYT, SILANET !jwrm2> unused, STOCKAAT1361: &                        NPYSITE, MULTISITEPYT, SILANET !jwrm2> unused, STOCKAAT
1393:       USE pymodule, ONLY : SITECOORDS,ELLST1,ELLMAT1362:       USE pymodule, ONLY : SITECOORDS,ELLST1,ELLMAT
1394: 1363: 
1395:       IMPLICIT NONE1364:       IMPLICIT NONE
1396:       1365:       
1397:       INTEGER :: J1, J2, J31366:       INTEGER :: J1, J2, J3
1398:       DOUBLE PRECISION :: XS(3*NTSITES)1367:       DOUBLE PRECISION :: XS(3*NTSITES)
1399:       CHARACTER(LEN=132),INTENT(IN) :: C1368:       CHARACTER(LEN=132),INTENT(IN) :: C
1400:       DOUBLE PRECISION :: Q(3*NATOMS),RMATBEST(3,3)1369:       DOUBLE PRECISION :: Q(3*NATOMS),RMATBEST(3,3)
1401: 1370: 
1402:       IF (DBPTDT ) THEN1371:       IF (DBPTDT ) THEN
1538:                     DO J3=1,NPYSITE1507:                     DO J3=1,NPYSITE
1539:                         WRITE(55,'(a5,2x,3f20.10,2x,a8,12f15.8,2x,a11,3f15.8)') 'O',&1508:                         WRITE(55,'(a5,2x,3f20.10,2x,a8,12f15.8,2x,a11,3f15.8)') 'O',&
1540:                         SITECOORDS(J3,1),SITECOORDS(J3,2),SITECOORDS(J3,3),&1509:                         SITECOORDS(J3,1),SITECOORDS(J3,2),SITECOORDS(J3,3),&
1541:                         'ellipse ',ELLST1(J3,1)*2.0D0,ELLST1(J3,2)*2.0D0,ELLST1(J3,3)*2.0D0,&1510:                         'ellipse ',ELLST1(J3,1)*2.0D0,ELLST1(J3,2)*2.0D0,ELLST1(J3,3)*2.0D0,&
1542:                         ELLMAT(J3,1,1),ELLMAT(J3,1,2),ELLMAT(J3,1,3),ELLMAT(J3,2,1),ELLMAT(J3,2,2),ELLMAT(J3,2,3),&1511:                         ELLMAT(J3,1,1),ELLMAT(J3,1,2),ELLMAT(J3,1,3),ELLMAT(J3,2,1),ELLMAT(J3,2,2),ELLMAT(J3,2,3),&
1543:                         ELLMAT(J3,3,1),ELLMAT(J3,3,2),ELLMAT(J3,3,3),&1512:                         ELLMAT(J3,3,1),ELLMAT(J3,3,2),ELLMAT(J3,3,3),&
1544:                         'atom_vector',Q(3*NATOMS/2+3*(J2-1)+1),&1513:                         'atom_vector',Q(3*NATOMS/2+3*(J2-1)+1),&
1545:                         Q(3*NATOMS/2+3*(J2-1)+2),Q(3*NATOMS/2+3*(J2-1)+3)1514:                         Q(3*NATOMS/2+3*(J2-1)+2),Q(3*NATOMS/2+3*(J2-1)+3)
1546:                     END DO1515:                     END DO
1547:              END DO1516:              END DO
1548:       ELSEIF (PYT) THEN ! sf344> implementation still buggy, use MULTISITEPY instead! 
1549:          WRITE(55,'(i6)') NATOMS*NRBSITES/2 
1550:          WRITE(55,'(1x,a)') trim(adjustl(c)) 
1551:         WRITE(*,*) 'in here' 
1552:          CALL PY_OUTPUT(55,Q) 
1553: !      ELSEIF (STOCKAAT) THEN1517: !      ELSEIF (STOCKAAT) THEN
1554: 1518: 
1555: !         WRITE(55,'(i6)') (NATOMS/2)*NRBSITES1519: !         WRITE(55,'(i6)') (NATOMS/2)*NRBSITES
1556: !         WRITE(55,'(1x,a)') trim(adjustl(c))1520: !         WRITE(55,'(1x,a)') trim(adjustl(c))
1557: !         DO J1 = 1, NATOMS/21521: !         DO J1 = 1, NATOMS/2
1558: !            J2 = (J1-1)*31522: !            J2 = (J1-1)*3
1559: !            WRITE(55,'(A5,1X,3F20.10)') 'O ', XS(J2+1), XS(J2+2), XS(J2+3)1523: !            WRITE(55,'(A5,1X,3F20.10)') 'O ', XS(J2+1), XS(J2+2), XS(J2+3)
1560: !         ENDDO1524: !         ENDDO
1561: 1525: 
1562:       ELSEIF (SILANET) THEN1526:       ELSEIF (SILANET) THEN
1600: 1564: 
1601:       ENDDO1565:       ENDDO
1602: 1566: 
1603:       END SUBROUTINE STFRAME1567:       END SUBROUTINE STFRAME
1604: 1568: 
1605: !     ----------------------------------------------------------------------------------------------1569: !     ----------------------------------------------------------------------------------------------
1606: 1570: 
1607:       SUBROUTINE RBROTMAT(T, X, ROTMAT, NTSITES)1571:       SUBROUTINE RBROTMAT(T, X, ROTMAT, NTSITES)
1608: !     takes the set of rigid-body coordinates T and returns X after rotation via the quaternion Q21572: !     takes the set of rigid-body coordinates T and returns X after rotation via the quaternion Q2
1609: !     about the origin1573: !     about the origin
1610:       USE KEY, ONLY : EFIELDT, MULTISITEPYT1574: 
1611:       IMPLICIT NONE1575:       IMPLICIT NONE
1612: 1576: 
1613:       INTEGER          :: I, J, NTSITES1577:       INTEGER          :: I, J, NTSITES
1614:       DOUBLE PRECISION :: T(3*NTSITES), X(3*NTSITES), T1(1:3), ROTMAT(3,3) !jwmr2> unused, Q2(4)1578:       DOUBLE PRECISION :: T(3*NTSITES), X(3*NTSITES), T1(1:3), ROTMAT(3,3) !jwmr2> unused, Q2(4)
1615:       DOUBLE PRECISION :: CMX, CMY, CMZ1579:       DOUBLE PRECISION :: CMX, CMY, CMZ
1616: !1580: !
1617: !     Move centre of mass to the origin.1581: !     Move centre of mass to the origin.
1618: !1582: !
1619:       CMX = 0.D0; CMY = 0.D0; CMZ = 0.D01583:       CMX = 0.D0; CMY = 0.D0; CMZ = 0.D0
1620:       DO I = 1, NTSITES1584:       DO I = 1, NTSITES
1621:          J = 3*I1585:          J = 3*I
1622:          CMX = CMX + T(J-2)1586:          CMX = CMX + T(J-2)
1623:          CMY = CMY + T(J-1)1587:          CMY = CMY + T(J-1)
1624:          CMZ = CMZ + T(J)1588:          CMZ = CMZ + T(J)
1625:       ENDDO1589:       ENDDO
1626:       CMX = CMX/NTSITES; CMY = CMY/NTSITES; CMZ = CMZ/NTSITES1590:       CMX = CMX/NTSITES; CMY = CMY/NTSITES; CMZ = CMZ/NTSITES
1627:       DO I = 1, NTSITES1591:       DO I = 1, NTSITES
1628:          J      = 3*I1592:          J      = 3*I
1629:          X(J-2) = T(J-2) - CMX1593:          X(J-2) = T(J-2) - CMX
1630:          X(J-1) = T(J-1) - CMY1594:          X(J-1) = T(J-1) - CMY
1631: ! do not move in the z direction for gravity field (currently implemented in MULTISITEPY)1595:          X(J)   = T(J)   - CMZ
1632:         IF(EFIELDT.AND.MULTISITEPYT) CMZ=0.0D0  
1633:         X(J)   = T(J)   - CMZ 
1634:       ENDDO1596:       ENDDO
1635: 1597: 
1636: !     extract the rotation matrix corresponding to rotation via Q1598: !     extract the rotation matrix corresponding to rotation via Q
1637: 1599: 
1638: !      CALL QROTMAT(Q2,RM)1600: !      CALL QROTMAT(Q2,RM)
1639: 1601: 
1640:       DO I = 1, NTSITES1602:       DO I = 1, NTSITES
1641: 1603: 
1642:          J        = 3*I1604:          J        = 3*I
1643:          T1(1:3)  = MATMUL(ROTMAT(:,:), X(J-2:J))1605:          T1(1:3)  = MATMUL(ROTMAT(:,:), X(J-2:J))


r31852/rigidb.f90 2017-01-30 12:30:12.736067298 +0000 r31851/rigidb.f90 2017-01-30 12:30:17.072126262 +0000
501:       USE KEY501:       USE KEY
502:       USE MODHESS502:       USE MODHESS
503: 503: 
504:       IMPLICIT NONE504:       IMPLICIT NONE
505: 505: 
506:       INTEGER            :: NATOMS, I, J, J1, J2506:       INTEGER            :: NATOMS, I, J, J1, J2
507:       DOUBLE PRECISION   :: Q(3*NATOMS), EV(3*NATOMS,6+NATOMS/2), NRMFCT(6+NATOMS/2)507:       DOUBLE PRECISION   :: Q(3*NATOMS), EV(3*NATOMS,6+NATOMS/2), NRMFCT(6+NATOMS/2)
508:       DOUBLE PRECISION   :: CMX, CMY, CMZ, THETA, THETA2, THETAH, DUMMY508:       DOUBLE PRECISION   :: CMX, CMY, CMZ, THETA, THETA2, THETAH, DUMMY
509: 509: 
510: !     INITIALIZE510: !     INITIALIZE
 511: 
511:       EV(:,:)   = 0.D0512:       EV(:,:)   = 0.D0
512:       NRMFCT(:) = 0.D0513:       NRMFCT(:) = 0.D0
513:       CMX       = 0.D0514:       CMX       = 0.D0
514:       CMY       = 0.D0515:       CMY       = 0.D0
515:       CMZ       = 0.D0516:       CMZ       = 0.D0
516: 517: 
517:       SHIFTED = .TRUE.518:       SHIFTED = .TRUE.
518:       IF (EFIELDT) THEN519:       IF (EFIELDT) THEN
519:           NZERO = 4520:           NZERO = 4
520:       ELSE521:       ELSE
607:               EV(J,6)    = - Q(J-2) + CMX608:               EV(J,6)    = - Q(J-2) + CMX
608:               EV(J1-2,6) = 0.5D0*Q(J1) + Q(J1-1)*Q(J1-2)/THETA2 - 0.5D0*Q(J1-1)*Q(J1-2)/(THETA*TAN(THETAH))609:               EV(J1-2,6) = 0.5D0*Q(J1) + Q(J1-1)*Q(J1-2)/THETA2 - 0.5D0*Q(J1-1)*Q(J1-2)/(THETA*TAN(THETAH))
609:               EV(J1-1,6) = THETAH/TAN(THETAH) + Q(J1-1)**2/THETA2 - 0.5D0*Q(J1-1)**2/(THETA*TAN(THETAH))610:               EV(J1-1,6) = THETAH/TAN(THETAH) + Q(J1-1)**2/THETA2 - 0.5D0*Q(J1-1)**2/(THETA*TAN(THETAH))
610:               EV(J1,6)   = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH))611:               EV(J1,6)   = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH))
611:               NRMFCT(6)  = NRMFCT(6) + EV(J-2,6)**2 + EV(J,6)**2 + EV(J1-2,6)**2 + EV(J1-1,6)**2 + EV(J1,6)**2612:               NRMFCT(6)  = NRMFCT(6) + EV(J-2,6)**2 + EV(J,6)**2 + EV(J1-2,6)**2 + EV(J1-1,6)**2 + EV(J1,6)**2
612: 613: 
613:             ENDIF614:             ENDIF
614: 615: 
615:          ENDIF616:          ENDIF
616: 617: 
617:          IF (GBT.OR.GBDT.OR.(PYGT.AND.UNIAXT).OR.STOCKAAT.OR.((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT)) THEN618:          IF (GBT.OR.GBDT.OR.(PYGT.AND.UNIAXT).OR.STOCKAAT.OR.((PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT)) THEN
618: 619: 
619: !     ROTATION ABOUT THE SYMMETRY AXIS620: !     ROTATION ABOUT THE SYMMETRY AXIS
620:             IF (THETA == 0.D0) THEN621:             IF (THETA == 0.D0) THEN
621: 622: 
622:                EV(J1,NZERO+I)  = 1.D0 623:                EV(J1,NZERO+I)  = 1.D0 
623:                NRMFCT(NZERO+I) = NRMFCT(NZERO+I) + EV(J1,NZERO+I)**2624:                NRMFCT(NZERO+I) = NRMFCT(NZERO+I) + EV(J1,NZERO+I)**2
624:             625:             
625:             ELSE 626:             ELSE 
626:                EV(J1-2,NZERO+I) = 0.5D0*Q(J1-1) + Q(J1-2)*Q(J1)/THETA2 - 0.5D0*Q(J1-2)*Q(J1)/(THETA*TAN(THETAH))627:                EV(J1-2,NZERO+I) = 0.5D0*Q(J1-1) + Q(J1-2)*Q(J1)/THETA2 - 0.5D0*Q(J1-2)*Q(J1)/(THETA*TAN(THETAH))
627:                EV(J1-1,NZERO+I) = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH))628:                EV(J1-1,NZERO+I) = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH))
628:                EV(J1,NZERO+I)   = THETAH*SIN(THETA) + Q(J1)*Q(J1)/THETA2 &629:                EV(J1,NZERO+I)   = THETAH*SIN(THETA) + Q(J1)*Q(J1)/THETA2 &
629:                                 + 0.5D0*(THETA*COS(THETA)-Q(J1)*Q(J1)/THETA)/TAN(THETAH)630:                                 + 0.5D0*(THETA*COS(THETA)-Q(J1)*Q(J1)/THETA)/TAN(THETAH)
630:                NRMFCT(NZERO+I)  = NRMFCT(NZERO+I) + EV(J1-2,NZERO+I)**2 + EV(J1-1,NZERO+I)**2 + EV(J1,NZERO+I)**2  631:                NRMFCT(NZERO+I)  = NRMFCT(NZERO+I) + EV(J1-2,NZERO+I)**2 + EV(J1-1,NZERO+I)**2 + EV(J1,NZERO+I)**2  
631:             ENDIF632:             ENDIF
632:        633:        
633:          ENDIF634:          ENDIF
634: 635: 
635:       ENDDO636:       ENDDO
636: 637: 
637:       IF (GBT.OR.GBDT.OR.(PYGT.AND.UNIAXT).OR.STOCKAAT.OR.((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT)) THEN638:       IF (GBT.OR.GBDT.OR.(PYGT.AND.UNIAXT).OR.STOCKAAT.OR.((PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT)) THEN
638:          NZERO = NZERO + NATOMS/2639:          NZERO = NZERO + NATOMS/2
639:       ELSE640:       ELSE
640:          NZERO = NZERO641:          NZERO = NZERO
641:       ENDIF642:       ENDIF
642: 643: 
643:       DO J = 1, NZERO644:       DO J = 1, NZERO
644: 645: 
645:          NRMFCT(J) = DSQRT(NRMFCT(J))646:          NRMFCT(J) = DSQRT(NRMFCT(J))
646:          EV(:,J)   = EV(:,J)/NRMFCT(J)647:          EV(:,J)   = EV(:,J)/NRMFCT(J)
647: 648: 
794:                EV(J,6)    = - Q(J-2) + CMX795:                EV(J,6)    = - Q(J-2) + CMX
795:                EV(J1-2,6) = 0.5D0*Q(J1) + Q(J1-1)*Q(J1-2)/THETA2 - 0.5D0*Q(J1-1)*Q(J1-2)/(THETA*TAN(THETAH))796:                EV(J1-2,6) = 0.5D0*Q(J1) + Q(J1-1)*Q(J1-2)/THETA2 - 0.5D0*Q(J1-1)*Q(J1-2)/(THETA*TAN(THETAH))
796:                EV(J1-1,6) = THETAH/TAN(THETAH) + Q(J1-1)**2/THETA2 - 0.5D0*Q(J1-1)**2/(THETA*TAN(THETAH))797:                EV(J1-1,6) = THETAH/TAN(THETAH) + Q(J1-1)**2/THETA2 - 0.5D0*Q(J1-1)**2/(THETA*TAN(THETAH))
797:                EV(J1,6)   = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH))798:                EV(J1,6)   = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH))
798:                NRMFCT(6)  = NRMFCT(6) + EV(J-2,6)**2 + EV(J,6)**2 + EV(J1-2,6)**2 + EV(J1-1,6)**2 + EV(J1,6)**2799:                NRMFCT(6)  = NRMFCT(6) + EV(J-2,6)**2 + EV(J,6)**2 + EV(J1-2,6)**2 + EV(J1-1,6)**2 + EV(J1,6)**2
799: 800: 
800:              ENDIF801:              ENDIF
801: 802: 
802:          ENDIF803:          ENDIF
803: 804: 
804:          IF (GBT.OR.GBDT.OR.(PYGT.AND.UNIAXT).OR.STOCKAAT.OR.((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT)) THEN805:          IF (GBT.OR.GBDT.OR.(PYGT.AND.UNIAXT).OR.STOCKAAT.OR.((PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT)) THEN
805: 806: 
806: !     ROTATION ABOUT THE SYMMETRY AXIS807: !     ROTATION ABOUT THE SYMMETRY AXIS
807:             IF (THETA == 0.D0) THEN808:             IF (THETA == 0.D0) THEN
808: 809: 
809:                EV(J1,NZERO+I)  = 1.D0810:                EV(J1,NZERO+I)  = 1.D0
810:                NRMFCT(NZERO+I) = NRMFCT(NZERO+I) + EV(J1,NZERO+I)**2811:                NRMFCT(NZERO+I) = NRMFCT(NZERO+I) + EV(J1,NZERO+I)**2
811: 812: 
812:             ELSE813:             ELSE
813:                EV(J1-2,NZERO+I) = 0.5D0*Q(J1-1) + Q(J1-2)*Q(J1)/THETA2 - 0.5D0*Q(J1-2)*Q(J1)/(THETA*TAN(THETAH))814:                EV(J1-2,NZERO+I) = 0.5D0*Q(J1-1) + Q(J1-2)*Q(J1)/THETA2 - 0.5D0*Q(J1-2)*Q(J1)/(THETA*TAN(THETAH))
814:                EV(J1-1,NZERO+I) = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH))815:                EV(J1-1,NZERO+I) = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH))
815:                EV(J1,NZERO+I)   = THETAH*SIN(THETA) + Q(J1)*Q(J1)/THETA2 &816:                EV(J1,NZERO+I)   = THETAH*SIN(THETA) + Q(J1)*Q(J1)/THETA2 &
816:                                 + 0.5D0*(THETA*COS(THETA)-Q(J1)*Q(J1)/THETA)/TAN(THETAH)817:                                 + 0.5D0*(THETA*COS(THETA)-Q(J1)*Q(J1)/THETA)/TAN(THETAH)
817:                NRMFCT(NZERO+I)  = NRMFCT(NZERO+I) + EV(J1-2,NZERO+I)**2 + EV(J1-1,NZERO+I)**2 + EV(J1,NZERO+I)**2818:                NRMFCT(NZERO+I)  = NRMFCT(NZERO+I) + EV(J1-2,NZERO+I)**2 + EV(J1-1,NZERO+I)**2 + EV(J1,NZERO+I)**2
818: 819: 
819:             ENDIF820:             ENDIF
820: 821: 
821:          ENDIF822:          ENDIF
822: 823: 
823:       ENDDO824:       ENDDO
824: 825: 
825:       IF (GBT.OR.GBDT.OR.(PYGT.AND.UNIAXT).OR.STOCKAAT.OR.((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT)) THEN826:       IF (GBT.OR.GBDT.OR.(PYGT.AND.UNIAXT).OR.STOCKAAT.OR.((PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT)) THEN
826:          NZERO = NZERO + NATOMS/2827:          NZERO = NZERO + NATOMS/2
827:       ELSE828:       ELSE
828:          NZERO = NZERO829:          NZERO = NZERO
829:       ENDIF830:       ENDIF
830: 831: 
831:       DO J = 1, NZERO  832:       DO J = 1, NZERO  
832: 833: 
833:          NRMFCT(J) = DSQRT(NRMFCT(J))834:          NRMFCT(J) = DSQRT(NRMFCT(J))
834:          EV(:,J)   = EV(:,J)/NRMFCT(J)835:          EV(:,J)   = EV(:,J)/NRMFCT(J)
835: 836: 
1163:       LOGICAL, INTENT(IN)          :: DEBUG1164:       LOGICAL, INTENT(IN)          :: DEBUG
1164:       DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)1165:       DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)
1165:       DOUBLE PRECISION dr, boxl(3)1166:       DOUBLE PRECISION dr, boxl(3)
1166:       INTEGER J1, J2, NSIZE1167:       INTEGER J1, J2, NSIZE
1167: 1168: 
1168:       !write(*,*) "RBMINDIST_BULK>"1169:       !write(*,*) "RBMINDIST_BULK>"
1169: 1170: 
1170:       BOXL(1) = PARAM11171:       BOXL(1) = PARAM1
1171:       BOXL(2) = PARAM21172:       BOXL(2) = PARAM2
1172:       BOXL(3) = PARAM31173:       BOXL(3) = PARAM3
1173: ! sf344 modification1174: 
1174: !      BOXL(:) = 100.0D0 
1175:       NSIZE = NTSITES1175:       NSIZE = NTSITES
1176:       ALLOCATE(XA(3*NSIZE),XB(3*NSIZE))1176:       ALLOCATE(XA(3*NSIZE),XB(3*NSIZE))
1177: 1177: 
1178: !     CONVERT TO THE CARTESIAN CORDINATES FOR THE SITES1178: !     CONVERT TO THE CARTESIAN CORDINATES FOR THE SITES
1179: 1179: 
1180:       CALL SITEPOS(RA,XA)1180:       CALL SITEPOS(RA,XA)
1181: 1181: 
1182:       CALL SITEPOS(RB,XB)1182:       CALL SITEPOS(RB,XB)
1183: 1183: 
1184: !     apply periodic boundary conditions1184: !     apply periodic boundary conditions
1210: 1210: 
1211:       SUBROUTINE RBMINDIST(RA,RB,NATOMS,DIST,Q2,DEBUG)1211:       SUBROUTINE RBMINDIST(RA,RB,NATOMS,DIST,Q2,DEBUG)
1212: 1212: 
1213: !     Follows the prescription of Kearsley, Acta Cryst. A, 45, 208-210, 1989, making necessary changes 1213: !     Follows the prescription of Kearsley, Acta Cryst. A, 45, 208-210, 1989, making necessary changes 
1214: !     to conform to right-handed rotation in the right-handed coordinate system.1214: !     to conform to right-handed rotation in the right-handed coordinate system.
1215: !     Brings RB to the best alignment with RA at the centre of mass of RA1215: !     Brings RB to the best alignment with RA at the centre of mass of RA
1216: !     Returns DIST as the actual distance, rather than the squared distance 1216: !     Returns DIST as the actual distance, rather than the squared distance 
1217: !     http://dx.doi.org/10.1107/S01087673880101281217: !     http://dx.doi.org/10.1107/S0108767388010128
1218: 1218: 
1219:       USE COMMONS, ONLY: NRBSITES1219:       USE COMMONS, ONLY: NRBSITES
1220:       USE KEY, ONLY: NTSITES, DBPT, DBPTDT, DMBLPYT, MSSTOCKT, STOCKAAT, EFIELDT, BULKT, MULTISITEPYT1220:       USE KEY, ONLY: NTSITES, DBPT, DBPTDT, DMBLPYT, MSSTOCKT, STOCKAAT, EFIELDT, BULKT
1221: 1221: 
1222:       IMPLICIT NONE1222:       IMPLICIT NONE
1223: 1223: 
1224:       INTEGER :: NATOMS1224:       INTEGER :: NATOMS
1225:       DOUBLE PRECISION :: RA(3*NATOMS)1225:       DOUBLE PRECISION :: RA(3*NATOMS)
1226:       DOUBLE PRECISION :: RB(3*NATOMS)1226:       DOUBLE PRECISION :: RB(3*NATOMS)
1227:       DOUBLE PRECISION :: DIST, Q2(4)1227:       DOUBLE PRECISION :: DIST, Q2(4)
1228:       LOGICAL          :: DEBUG1228:       LOGICAL          :: DEBUG
1229:       INTEGER          :: J1, J2, NSIZE, JMIN, INFO1229:       INTEGER          :: J1, J2, NSIZE, JMIN, INFO
1230:       DOUBLE PRECISION :: QMAT(4,4), TEMPA(9*NATOMS), XM, YM, ZM, XP, YP, ZP1230:       DOUBLE PRECISION :: QMAT(4,4), TEMPA(9*NATOMS), XM, YM, ZM, XP, YP, ZP
1231:       DOUBLE PRECISION :: DIAG(4), MINV, CMXA, CMYA, CMZA, CMXB, CMYB, CMZB1231:       DOUBLE PRECISION :: DIAG(4), MINV, CMXA, CMYA, CMZA, CMXB, CMYB, CMZB
1232:       DOUBLE PRECISION :: R(3), P(3), RM(3,3)1232:       DOUBLE PRECISION :: R(3), P(3), RM(3,3)
1233:       DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)1233:       DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)
1234:       DOUBLE PRECISION :: ENERGY, VNEW(3*NATOMS), RMS1234:       DOUBLE PRECISION :: ENERGY, VNEW(3*NATOMS), RMS
1235: 1235: 
1236:       IF ((DBPT.AND.EFIELDT).OR.(DBPTDT.AND.EFIELDT).OR.(DMBLPYT.AND.EFIELDT)  &1236:       IF ((DBPT.AND.EFIELDT).OR.(DBPTDT.AND.EFIELDT).OR.(DMBLPYT.AND.EFIELDT)  &
1237:      &    .OR.(MSSTOCKT.AND.EFIELDT).OR.(STOCKAAT.AND.EFIELDT).OR.(MULTISITEPYT.AND.EFIELDT)) THEN1237:      &    .OR.(MSSTOCKT.AND.EFIELDT).OR.(STOCKAAT.AND.EFIELDT)) THEN
1238: 1238: 
1239:          CALL FLDMINDIST(RA,RB,NATOMS,DIST,DEBUG,Q2)1239:          CALL FLDMINDIST(RA,RB,NATOMS,DIST,DEBUG,Q2)
1240:          RETURN1240:          RETURN
1241:     1241: 
1242:       ENDIF 1242:       ENDIF 
1243: 1243: 
1244:       IF (BULKT) THEN1244:       IF (BULKT) THEN
1245:          CALL RBMINDIST_BULK(RA,RB,NATOMS,DIST,Q2,DEBUG)1245:          CALL RBMINDIST_BULK(RA,RB,NATOMS,DIST,Q2,DEBUG)
1246:          RETURN1246:          RETURN
1247:       ENDIF1247:       ENDIF
1248: 1248: 
1249:       NSIZE = NTSITES1249:       NSIZE = NTSITES
1250:       ALLOCATE(XA(3*NSIZE),XB(3*NSIZE))1250:       ALLOCATE(XA(3*NSIZE),XB(3*NSIZE))
1251: 1251: 
1271:          J2 = 3*(J1-1)1271:          J2 = 3*(J1-1)
1272:          CMXB = CMXB + RB(J2+1)1272:          CMXB = CMXB + RB(J2+1)
1273:          CMYB = CMYB + RB(J2+2)1273:          CMYB = CMYB + RB(J2+2)
1274:          CMZB = CMZB + RB(J2+3)1274:          CMZB = CMZB + RB(J2+3)
1275:       ENDDO1275:       ENDDO
1276:       CMXB = 2.D0*CMXB/NATOMS; CMYB = 2.D0*CMYB/NATOMS; CMZB = 2.D0*CMZB/NATOMS1276:       CMXB = 2.D0*CMXB/NATOMS; CMYB = 2.D0*CMYB/NATOMS; CMZB = 2.D0*CMZB/NATOMS
1277:       DO J1 = 1, NATOMS/21277:       DO J1 = 1, NATOMS/2
1278:          J2 = 3*(J1-1)1278:          J2 = 3*(J1-1)
1279:          RB(J2+1) = RB(J2+1) - CMXB1279:          RB(J2+1) = RB(J2+1) - CMXB
1280:          RB(J2+2) = RB(J2+2) - CMYB1280:          RB(J2+2) = RB(J2+2) - CMYB
1281: !sf344 modification          
1282:          RB(J2+3) = RB(J2+3) - CMZB1281:          RB(J2+3) = RB(J2+3) - CMZB
1283:       ENDDO1282:       ENDDO
1284: 1283: 
1285: !     CONVERT TO THE CARTESIAN CORDINATES FOR THE SITES1284: !     CONVERT TO THE CARTESIAN CORDINATES FOR THE SITES
1286: 1285: 
1287:       CALL SITEPOS(RA,XA)1286:       CALL SITEPOS(RA,XA)
1288:       CALL SITEPOS(RB,XB)1287:       CALL SITEPOS(RB,XB)
1289:  1288:  
1290: !     Create matrix QMAT1289: !     Create matrix QMAT
1291: 1290: 
1340:       DIST = SQRT(MINV)1339:       DIST = SQRT(MINV)
1341: 1340: 
1342: !     Calculate Q2 from QMAT and JMIN1341: !     Calculate Q2 from QMAT and JMIN
1343:       Q2(1) = QMAT(1,JMIN); Q2(2) = QMAT(2,JMIN); Q2(3) = QMAT(3,JMIN); Q2(4) = QMAT(4,JMIN)1342:       Q2(1) = QMAT(1,JMIN); Q2(2) = QMAT(2,JMIN); Q2(3) = QMAT(3,JMIN); Q2(4) = QMAT(4,JMIN)
1344: 1343: 
1345: !     Update RB based on center of mass1344: !     Update RB based on center of mass
1346:       DO J1 = 1, NATOMS/21345:       DO J1 = 1, NATOMS/2
1347:          J2 = 3*(J1-1)1346:          J2 = 3*(J1-1)
1348:          RB(J2+1) = RB(J2+1) - CMXB1347:          RB(J2+1) = RB(J2+1) - CMXB
1349:          RB(J2+2) = RB(J2+2) - CMYB1348:          RB(J2+2) = RB(J2+2) - CMYB
1350: ! sf344 modification          
1351:          RB(J2+3) = RB(J2+3) - CMZB1349:          RB(J2+3) = RB(J2+3) - CMZB
1352:       ENDDO1350:       ENDDO
1353: 1351: 
1354: !     rotate RB based on Q2 and center of mass1352: !     rotate RB based on Q2 and center of mass
1355:       CALL RBNEWROTGEOM(NATOMS,RB,Q2,RM,CMXA,CMYA,CMZA)1353:       CALL RBNEWROTGEOM(NATOMS,RB,Q2,RM,CMXA,CMYA,CMZA)
1356: 1354: 
1357:       !write(*,*) "RBMINDIST> after  ", RA(1), RA(2), RA(3), RB(1), RB(2), RB(3)1355:       !write(*,*) "RBMINDIST> after  ", RA(1), RA(2), RA(3), RB(1), RB(2), RB(3)
1358:       DEALLOCATE(XA,XB)1356:       DEALLOCATE(XA,XB)
1359: 1357: 
1360:       END SUBROUTINE RBMINDIST1358:       END SUBROUTINE RBMINDIST
1376: 1374: 
1377:       CALL QROTMAT(Q2,RM)1375:       CALL QROTMAT(Q2,RM)
1378: 1376: 
1379:       DO I = 1, NATOMS/21377:       DO I = 1, NATOMS/2
1380:   1378:   
1381:          J    = 3*(I-1)1379:          J    = 3*(I-1)
1382:          R(:) = MATMUL(RM(:,:), COORDS(J+1:J+3))1380:          R(:) = MATMUL(RM(:,:), COORDS(J+1:J+3))
1383: 1381: 
1384:          COORDS(J+1) = R(1) + CX1382:          COORDS(J+1) = R(1) + CX
1385:          COORDS(J+2) = R(2) + CY1383:          COORDS(J+2) = R(2) + CY
1386:          COORDS(J+3) = R(3) !+ CZ1384:          COORDS(J+3) = R(3) + CZ
1387:       1385:       
1388: !     CONVERT THE ANGLE-AXIS COORDINATES1386: !     CONVERT THE ANGLE-AXIS COORDINATES
1389: 1387: 
1390:          J      = 3*NATOMS/2 + J1388:          J      = 3*NATOMS/2 + J
1391:          P(:)   = COORDS(J+1:J+3)1389:          P(:)   = COORDS(J+1:J+3)
1392: 1390: 
1393:          CALL QROTAA(Q2,P)1391:          CALL QROTAA(Q2,P)
1394: 1392: 
1395:          COORDS(J+1:J+3) = P(1:3)1393:          COORDS(J+1:J+3) = P(1:3)
1396: 1394: 
1398: 1396: 
1399:       END SUBROUTINE RBNEWROTGEOM1397:       END SUBROUTINE RBNEWROTGEOM
1400: 1398: 
1401: !     ----------------------------------------------------------------------------------------------1399: !     ----------------------------------------------------------------------------------------------
1402: 1400: 
1403:       SUBROUTINE FLDMINDIST(RA,RB,NATOMS,DIST,DEBUG,Q2)1401:       SUBROUTINE FLDMINDIST(RA,RB,NATOMS,DIST,DEBUG,Q2)
1404: 1402: 
1405: !     returns DIST as the actual distance, rather than the squared distance1403: !     returns DIST as the actual distance, rather than the squared distance
1406: 1404: 
1407:       USE COMMONS, ONLY: NRBSITES1405:       USE COMMONS, ONLY: NRBSITES
1408:       USE KEY, ONLY: NTSITES, STOCKAAT, MULTISITEPYT1406:       USE KEY, ONLY: NTSITES, STOCKAAT
1409: 1407: 
1410:       IMPLICIT NONE1408:       IMPLICIT NONE
1411: 1409: 
1412:       INTEGER          :: J1, J2, J3, J4, NATOMS, NSIZE, JMIN, INFO1410:       INTEGER          :: J1, J2, J3, J4, NATOMS, NSIZE, JMIN, INFO
1413:       DOUBLE PRECISION :: RA(3*NATOMS), RB(3*NATOMS), DIST, QMAT(2,2), XM, YM, ZM, XP, YP, ZP1411:       DOUBLE PRECISION :: RA(3*NATOMS), RB(3*NATOMS), DIST, QMAT(2,2), XM, YM, ZM, XP, YP, ZP
1414:       DOUBLE PRECISION :: MINV, Q2(4), CMXA, CMYA, CMZA, CMXB, CMYB, CMZB1412:       DOUBLE PRECISION :: MINV, Q2(4), CMXA, CMYA, CMZA, CMXB, CMYB, CMZB
1415:       DOUBLE PRECISION :: R(3), P(3), RM(3,3) 1413:       DOUBLE PRECISION :: R(3), P(3), RM(3,3) 
1416:       DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)1414:       DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)
1417:       DOUBLE PRECISION :: ENERGY, VNEW(3*NATOMS), RMS, DUMMY1415:       DOUBLE PRECISION :: ENERGY, VNEW(3*NATOMS), RMS, DUMMY
1418:       LOGICAL          :: DEBUG1416:       LOGICAL          :: DEBUG
1427:          J2 = 3*(J1-1)1425:          J2 = 3*(J1-1)
1428:          CMXA = CMXA + RA(J2+1)1426:          CMXA = CMXA + RA(J2+1)
1429:          CMYA = CMYA + RA(J2+2)1427:          CMYA = CMYA + RA(J2+2)
1430:          CMZA = CMZA + RA(J2+3)1428:          CMZA = CMZA + RA(J2+3)
1431:       ENDDO1429:       ENDDO
1432:       CMXA = 2.D0*CMXA/NATOMS; CMYA = 2.D0*CMYA/NATOMS; CMZA = 2.D0*CMZA/NATOMS1430:       CMXA = 2.D0*CMXA/NATOMS; CMYA = 2.D0*CMYA/NATOMS; CMZA = 2.D0*CMZA/NATOMS
1433:       DO J1 = 1, NATOMS/21431:       DO J1 = 1, NATOMS/2
1434:          J2 = 3*(J1-1)1432:          J2 = 3*(J1-1)
1435:          RA(J2+1) = RA(J2+1) - CMXA1433:          RA(J2+1) = RA(J2+1) - CMXA
1436:          RA(J2+2) = RA(J2+2) - CMYA1434:          RA(J2+2) = RA(J2+2) - CMYA
1437: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)1435:          RA(J2+3) = RA(J2+3) - CMZA
1438:          IF(MULTISITEPYT) CMZA=0.0D0  
1439:          RA(J2+3) = RA(J2+3) - CMZA  
1440:       ENDDO1436:       ENDDO
1441: 1437: 
1442:       CMXB = 0.0D0; CMYB = 0.0D0; CMZB = 0.0D01438:       CMXB = 0.0D0; CMYB = 0.0D0; CMZB = 0.0D0
1443:       DO J1 = 1, NATOMS/21439:       DO J1 = 1, NATOMS/2
1444:          J2 = 3*(J1-1)1440:          J2 = 3*(J1-1)
1445:          CMXB = CMXB + RB(J2+1)1441:          CMXB = CMXB + RB(J2+1)
1446:          CMYB = CMYB + RB(J2+2)1442:          CMYB = CMYB + RB(J2+2)
1447:          CMZB = CMZB + RB(J2+3)1443:          CMZB = CMZB + RB(J2+3)
1448:       ENDDO1444:       ENDDO
1449:       CMXB = 2.D0*CMXB/NATOMS; CMYB = 2.D0*CMYB/NATOMS; CMZB = 2.D0*CMZB/NATOMS1445:       CMXB = 2.D0*CMXB/NATOMS; CMYB = 2.D0*CMYB/NATOMS; CMZB = 2.D0*CMZB/NATOMS
1450:       DO J1 = 1, NATOMS/21446:       DO J1 = 1, NATOMS/2
1451:          J2 = 3*(J1-1)1447:          J2 = 3*(J1-1)
1452:          RB(J2+1) = RB(J2+1) - CMXB1448:          RB(J2+1) = RB(J2+1) - CMXB
1453:          RB(J2+2) = RB(J2+2) - CMYB1449:          RB(J2+2) = RB(J2+2) - CMYB
1454: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY) 
1455:          IF(MULTISITEPYT) CMZB=0.0D0 
1456:          RB(J2+3) = RB(J2+3) - CMZB1450:          RB(J2+3) = RB(J2+3) - CMZB
1457:       ENDDO1451:       ENDDO
1458: 1452: 
1459: !     CONVERT TO THE CARTESIAN CORDINATES FOR THE SITES1453: !     CONVERT TO THE CARTESIAN CORDINATES FOR THE SITES
1460: 1454: 
1461:       IF (STOCKAAT) THEN1455:       IF (STOCKAAT) THEN
1462:          XA(:) = RA(:)1456:          XA(:) = RA(:)
1463:          XB(:) = RB(:)1457:          XB(:) = RB(:)
1464:       ELSE1458:       ELSE
1465:          CALL SITEPOS(RA,XA)1459:          CALL SITEPOS(RA,XA)
1517:             PRINT '(A,G20.10,A)','newmindist> WARNING MINV is ',MINV,' change to absolute value'1511:             PRINT '(A,G20.10,A)','newmindist> WARNING MINV is ',MINV,' change to absolute value'
1518:             MINV = -MINV1512:             MINV = -MINV
1519:          ENDIF1513:          ENDIF
1520:       ENDIF1514:       ENDIF
1521:       DIST = SQRT(MINV)1515:       DIST = SQRT(MINV)
1522:   1516:   
1523:       DO J1 = 1, NATOMS/21517:       DO J1 = 1, NATOMS/2
1524:          J2 = 3*(J1-1)1518:          J2 = 3*(J1-1)
1525:          RB(J2+1) = RB(J2+1) - CMXB1519:          RB(J2+1) = RB(J2+1) - CMXB
1526:          RB(J2+2) = RB(J2+2) - CMYB1520:          RB(J2+2) = RB(J2+2) - CMYB
1527: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY) 
1528:          IF(MULTISITEPYT) CMZB=0.0D0  
1529:          RB(J2+3) = RB(J2+3) - CMZB1521:          RB(J2+3) = RB(J2+3) - CMZB
1530:       ENDDO1522:       ENDDO
1531: 1523: 
1532:       CALL RBNEWROTGEOM(NATOMS,RB,Q2,RM,CMXA,CMYA,CMZA)1524:       CALL RBNEWROTGEOM(NATOMS,RB,Q2,RM,CMXA,CMYA,CMZA)
1533: 1525: 
1534:       DEALLOCATE(XA,XB)1526:       DEALLOCATE(XA,XB)
1535: 1527: 
1536:       END SUBROUTINE FLDMINDIST1528:       END SUBROUTINE FLDMINDIST
1537: 1529: 
1538: !     ----------------------------------------------------------------------------------------------1530: !     ----------------------------------------------------------------------------------------------
1650: !     ----------------------------------------------------------------------------------------------1642: !     ----------------------------------------------------------------------------------------------
1651: 1643: 
1652: 1644: 
1653: SUBROUTINE UNIAXGETPATHLENGTH(RA,RB,TEMP)1645: SUBROUTINE UNIAXGETPATHLENGTH(RA,RB,TEMP)
1654: use key, only : PYA11646: use key, only : PYA1
1655: use commons, only : NATOMS1647: use commons, only : NATOMS
1656: implicit none1648: implicit none
1657: 1649: 
1658: integer          :: NSIZE, J1, J2, J3, J4, NRBSITES1650: integer          :: NSIZE, J1, J2, J3, J4, NRBSITES
1659: double precision :: XA(9*NATOMS/2), XB(9*NATOMS/2), RA(3*NATOMS), RB(3*NATOMS), P(3), RM(3,3), RBSITE(3,3), TEMP1651: double precision :: XA(9*NATOMS/2), XB(9*NATOMS/2), RA(3*NATOMS), RB(3*NATOMS), P(3), RM(3,3), RBSITE(3,3), TEMP
1660: double precision :: R(3),rbdistance1652: double precision :: R(3)
1661: 1653: 
1662:       NRBSITES = 31654:       NRBSITES = 3
1663:       NSIZE = NRBSITES*NATOMS/21655:       NSIZE = NRBSITES*NATOMS/2
1664: ! sf344> try to generalize this for ellipsoids, using sites displaced by 1 along the symmetry axis       
1665:       IF(PYA1(3)==0.0D0) THEN 
1666:         rbdistance= 1.0D0 
1667:       ELSE 
1668:         rbdistance=PYA1(3) 
1669:       END IF 
1670: 1656: 
1671:       RBSITE(1,1) = 0.0D01657:       RBSITE(1,1) = 0.0D0
1672:       RBSITE(1,2) = 0.0D01658:       RBSITE(1,2) = 0.0D0
1673:       RBSITE(1,3) = 0.0D01659:       RBSITE(1,3) = 0.0D0
1674:       RBSITE(2,1) = 0.0D01660:       RBSITE(2,1) = 0.0D0
1675:       RBSITE(2,2) = 0.0D01661:       RBSITE(2,2) = 0.0D0
1676:       RBSITE(2,3) = rbdistance1662:       RBSITE(2,3) = PYA1(3)
1677:       RBSITE(3,1) = 0.0D01663:       RBSITE(3,1) = 0.0D0
1678:       RBSITE(3,2) = 0.0D01664:       RBSITE(3,2) = 0.0D0
1679:       RBSITE(3,3) = -rbdistance1665:       RBSITE(3,3) =-PYA1(3)
 1666: 
1680: !      ALLOCATE(XA(3*NSIZE),XB(3*NSIZE))1667: !      ALLOCATE(XA(3*NSIZE),XB(3*NSIZE))
1681: 1668: 
1682: !      write(*,*) 'ra'1669: !      write(*,*) 'ra'
1683: !      write(*,*) ra(:)1670: !      write(*,*) ra(:)
1684: !     ----------------------------------------------------------------------------------------------1671: !     ----------------------------------------------------------------------------------------------
1685: !      CALL LWOTPGH(RB,VNEW,ENERGY,.TRUE.,.FALSE.)1672: !      CALL LWOTPGH(RB,VNEW,ENERGY,.TRUE.,.FALSE.)
1686: !      WRITE(*,*) ENERGY, SQRT(DOT_PRODUCT(VNEW(:),VNEW(:)))1673: !      WRITE(*,*) ENERGY, SQRT(DOT_PRODUCT(VNEW(:),VNEW(:)))
1687: !      CALL DUMBBELLP(RB,VNEW,ENERGY,.TRUE.,.FALSE.)1674: !      CALL DUMBBELLP(RB,VNEW,ENERGY,.TRUE.,.FALSE.)
1688: !      WRITE(*,*) ENERGY, SQRT(DOT_PRODUCT(VNEW(:),VNEW(:)))1675: !      WRITE(*,*) ENERGY, SQRT(DOT_PRODUCT(VNEW(:),VNEW(:)))
1689: !     ----------------------------------------------------------------------------------------------1676: !     ----------------------------------------------------------------------------------------------
2174: 2161: 
2175:       IF (EIGENVECTORT) THEN      2162:       IF (EIGENVECTORT) THEN      
2176:          HESS = AP2163:          HESS = AP
2177:       ENDIF2164:       ENDIF
2178: 2165: 
2179:       END SUBROUTINE NRMLMD 2166:       END SUBROUTINE NRMLMD 
2180: 2167: 
2181: 2168: 
2182:       SUBROUTINE COMPUTEINERTIA(RMI, KBLOCK, TMASS) 2169:       SUBROUTINE COMPUTEINERTIA(RMI, KBLOCK, TMASS) 
2183: 2170: 
2184:         USE KEY, ONLY: NTIPT, PAPT, MULTISITEPYT2171:         USE KEY, ONLY: NTIPT, PAPT
2185:         IMPLICIT NONE2172:         IMPLICIT NONE
2186:         DOUBLE PRECISION :: TMASS, KBLOCK(3,3), RMI(3,3)       2173:         DOUBLE PRECISION :: TMASS, KBLOCK(3,3), RMI(3,3)       
2187: 2174: 
2188:         IF (NTIPT) THEN2175:         IF (NTIPT) THEN
2189:            2176:            
2190:            CALL INERTIANTIP(RMI, KBLOCK, TMASS)2177:            CALL INERTIANTIP(RMI, KBLOCK, TMASS)
2191:            2178:            
2192:         ELSE IF (PAPT) THEN2179:         ELSE IF (PAPT) THEN
2193: !jwrm2> added PAP potential2180: !jwrm2> added PAP potential
2194:            CALL INERTIAPAP(RMI, KBLOCK, TMASS)2181:            CALL INERTIAPAP(RMI, KBLOCK, TMASS)
2195:         ELSE IF (MULTISITEPYT) THEN2182: 
2196: !sf344> added PY potential (not generic yet) 
2197:            CALL INERTIAPY(RMI, KBLOCK, TMASS) 
2198:         
2199:         ENDIF2183:         ENDIF
2200: 2184: 
2201:       END SUBROUTINE COMPUTEINERTIA2185:       END SUBROUTINE COMPUTEINERTIA


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0