hdiff output

r32675/potential.f90 2017-06-01 17:30:08.429289121 +0100 r32674/potential.f90 2017-06-01 17:30:08.649296507 +0100
907:          CSMGPSAVE=CSMGP907:          CSMGPSAVE=CSMGP
908:          CSMGP=CSMGUIDEGP908:          CSMGP=CSMGUIDEGP
909:          CSMGPINDEXSAVE=CSMGPINDEX909:          CSMGPINDEXSAVE=CSMGPINDEX
910:          CSMGPINDEX=CSMGUIDEGPINDEX910:          CSMGPINDEX=CSMGUIDEGPINDEX
911:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPGUIDE(1:3, 1:3, 1:2*CSMGPINDEX)911:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPGUIDE(1:3, 1:3, 1:2*CSMGPINDEX)
912:       END IF912:       END IF
913:       IF (PERMDIST) THEN913:       IF (PERMDIST) THEN
914:          DO J1=1, CSMGPINDEX914:          DO J1=1, CSMGPINDEX
915:             XTEMP(1:3*NATOMS)=X(1:3*NATOMS)915:             XTEMP(1:3*NATOMS)=X(1:3*NATOMS)
916:             CALL CSMROT(XTEMP, DUMMY, 1, J1)916:             CALL CSMROT(XTEMP, DUMMY, 1, J1)
917:             CALL ALIGN_DECIDE(XTEMP, DUMMY, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)917:             CALL MINPERMDIST(XTEMP, DUMMY, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)
918:             CALL CSMROT(DUMMY, XTEMP,-1, J1) ! need to rotate the permuted rotated images back to the reference orientation918:             CALL CSMROT(DUMMY, XTEMP,-1, J1) ! need to rotate the permuted rotated images back to the reference orientation
919:             CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)919:             CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)
920:          END DO920:          END DO
921:       ELSE921:       ELSE
922:          DO J1=1, CSMGPINDEX922:          DO J1=1, CSMGPINDEX
923:             CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=X(1:3*NATOMS)923:             CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=X(1:3*NATOMS)
924:          END DO924:          END DO
925:       END IF925:       END IF
926:       CALL CSMMIN(X, EREAL, CSMRMS, CSMIT)926:       CALL CSMMIN(X, EREAL, CSMRMS, CSMIT)
927: 927: 
958:          END DO958:          END DO
959:          CSMAV(1:3*NATOMS)=CSMAV(1:3*NATOMS)/CSMGPINDEX959:          CSMAV(1:3*NATOMS)=CSMAV(1:3*NATOMS)/CSMGPINDEX
960: !960: !
961: !  Check the CSM for the averaged structure. It should be zero if this structure has the961: !  Check the CSM for the averaged structure. It should be zero if this structure has the
962: !  right point group. Need to reset CSMIMAGES and CSMNORM temporarily.962: !  right point group. Need to reset CSMIMAGES and CSMNORM temporarily.
963: !963: !
964:          IF (PERMDIST) THEN964:          IF (PERMDIST) THEN
965:             DO J1=1, CSMGPINDEX965:             DO J1=1, CSMGPINDEX
966:                XTEMP(1:3*NATOMS)=CSMAV(1:3*NATOMS)966:                XTEMP(1:3*NATOMS)=CSMAV(1:3*NATOMS)
967:                CALL CSMROT(XTEMP, DUMMY, 1, J1)967:                CALL CSMROT(XTEMP, DUMMY, 1, J1)
968:                CALL ALIGN_DECIDE(XTEMP, DUMMY, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, DUMMY2, DIST2, RIGID, RMAT)968:                CALL MINPERMDIST(XTEMP, DUMMY, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, DUMMY2, DIST2, RIGID, RMAT)
969:                CALL CSMROT(DUMMY, XTEMP,-1, J1) ! need to rotate the permuted rotated images back to the reference orientation969:                CALL CSMROT(DUMMY, XTEMP,-1, J1) ! need to rotate the permuted rotated images back to the reference orientation
970:                CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)970:                CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)
971:             END DO971:             END DO
972:          ELSE972:          ELSE
973:             DO J1=1, CSMGPINDEX973:             DO J1=1, CSMGPINDEX
974:                CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=CSMAV(1:3*NATOMS)974:                CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=CSMAV(1:3*NATOMS)
975:             END DO975:             END DO
976:          END IF976:          END IF
977:          AA(1)=0.0D0; AA(2)=0.0D0; AA(3)=6.283185307D0 ! should give an identity matrix977:          AA(1)=0.0D0; AA(2)=0.0D0; AA(3)=6.283185307D0 ! should give an identity matrix
978:          SAVECSMPMAT(1:3, 1:3)=CSMPMAT(1:3, 1:3)978:          SAVECSMPMAT(1:3, 1:3)=CSMPMAT(1:3, 1:3)
990:       END IF ! DEBUG990:       END IF ! DEBUG
991:       IF (CSMDOGUIDET) THEN ! undo guiding changes991:       IF (CSMDOGUIDET) THEN ! undo guiding changes
992:          CSMGP=CSMGPSAVE992:          CSMGP=CSMGPSAVE
993:          CSMGPINDEX=CSMGPINDEXSAVE993:          CSMGPINDEX=CSMGPINDEXSAVE
994:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPSAVE(1:3, 1:3, 1:2*CSMGPINDEX)994:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPSAVE(1:3, 1:3, 1:2*CSMGPINDEX)
995:          CSMNORM=CSMNORMSAVE995:          CSMNORM=CSMNORMSAVE
996:       END IF996:       END IF
997: 997: 
998:    ELSE IF (PERMOPT .OR. PERMINVOPT .OR. DISTOPT) THEN998:    ELSE IF (PERMOPT .OR. PERMINVOPT .OR. DISTOPT) THEN
999: !  EREAL is the distance in this case999: !  EREAL is the distance in this case
1000:       CALL ALIGN_DECIDE(FIN, X, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)1000:       CALL MINPERMDIST(FIN, X, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)
1001: 1001: 
1002:    ELSE IF (MKTRAPT) THEN1002:    ELSE IF (MKTRAPT) THEN
1003:       CALL MKTRAP(NATOMS, X, GRAD, EREAL)1003:       CALL MKTRAP(NATOMS, X, GRAD, EREAL)
1004:    ELSE IF (BLJCLUSTER) THEN1004:    ELSE IF (BLJCLUSTER) THEN
1005:       CALL RAD(X, GRAD, EREAL, GRADT)1005:       CALL RAD(X, GRAD, EREAL, GRADT)
1006:       CALL LJPSHIFTBINC(X, GRAD, EREAL, GRADT, SECT)1006:       CALL LJPSHIFTBINC(X, GRAD, EREAL, GRADT, SECT)
1007: 1007: 
1008:    ELSE IF (BLJCLUSTER_NOCUT) THEN1008:    ELSE IF (BLJCLUSTER_NOCUT) THEN
1009: ! ds656> Binary LJ without cutoff1009: ! ds656> Binary LJ without cutoff
1010:       CALL RAD(X, GRAD, EREAL, GRADT)1010:       CALL RAD(X, GRAD, EREAL, GRADT)
1345:       STOP1345:       STOP
1346: !      CALL ENERGY_3D_APBC(X, GRAD, EREAL, GRADT, SECT)1346: !      CALL ENERGY_3D_APBC(X, GRAD, EREAL, GRADT, SECT)
1347:    ELSE IF (THREEDPBCT) THEN1347:    ELSE IF (THREEDPBCT) THEN
1348:       WRITE(*, *) 'This potential option is commented out.'1348:       WRITE(*, *) 'This potential option is commented out.'
1349:       STOP1349:       STOP
1350: !      CALL ENERGY_3D_PBC(X, GRAD, EREAL, GRADT, SECT)1350: !      CALL ENERGY_3D_PBC(X, GRAD, EREAL, GRADT, SECT)
1351: 1351: 
1352:    ELSE1352:    ELSE
1353: ! Otherwise, assume Lennard-Jones type particles.1353: ! Otherwise, assume Lennard-Jones type particles.
1354: 1354: 
1355:       IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN 
1356:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS) 
1357:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS) 
1358:       END IF 
1359:  
1360: ! RAD must be called before the routine that calculates the potential or LBFGS1355: ! RAD must be called before the routine that calculates the potential or LBFGS
1361: ! will get confused even if EVAP is set .TRUE. correctly.1356: ! will get confused even if EVAP is set .TRUE. correctly.
1362:       CALL RAD(X, GRAD, EREAL, GRADT)1357:       CALL RAD(X, GRAD, EREAL, GRADT)
1363:       IF (EVAPREJECT) RETURN1358:       IF (EVAPREJECT) RETURN
1364:       IF (CUTT) THEN1359:       IF (CUTT) THEN
1365:          CALL LJCUT(X, GRAD, EREAL, GRADT, SECT)1360:          CALL LJCUT(X, GRAD, EREAL, GRADT, SECT)
1366:       ELSE1361:       ELSE
1367:          IF (CUDAT) THEN1362:          IF (CUDAT) THEN
1368: ! This call copies CPU coordinates to GPU, calculates energy/gradient and copies energy/gradient back to CPU1363: ! This call copies CPU coordinates to GPU, calculates energy/gradient and copies energy/gradient back to CPU
1369:             CALL CUDA_ENEGRAD_WRAPPER(NATOMS, X, EREAL, GRADATOMS)1364:             CALL CUDA_ENEGRAD_WRAPPER(NATOMS, X, EREAL, GRADATOMS)
1370:             GRAD(1:3*NATOMS)=GRADATOMS(:)1365:             GRAD(1:3*NATOMS)=GRADATOMS(:)
1371:          ELSE1366:          ELSE
1372:             CALL LJ(X, GRAD, EREAL, GRADT, SECT)1367:             CALL LJ(X, GRAD, EREAL, GRADT, SECT)
1373:          END IF1368:          END IF
1374:          IF (AXTELL) CALL AXT(NATOMS, X, GRAD, EREAL, GRADT, ZSTAR)1369:          IF (AXTELL) CALL AXT(NATOMS, X, GRAD, EREAL, GRADT, ZSTAR)
1375:       END IF1370:       END IF
1376:  
1377:        IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN 
1378:  
1379:            X(DEGFREEDOMS+1:3*NATOMS)=0.0D0 
1380:            X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS) 
1381:  
1382:            IF(GRADT) THEN 
1383:               GRADATOMS(:) = GRAD(:)  ! This is needed for AACONVERGENCET 
1384:               CALL TRANSFORMGRAD(GRAD, XRIGIDCOORDS, XRIGIDGRAD) 
1385:               GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0 
1386:               GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS) 
1387:            ENDIF 
1388:  
1389:            IF(SECT) THEN 
1390:               CALL TRANSFORMHESSIAN(HESS, GRAD, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET) 
1391:               HESS(DEGFREEDOMS+1:3*NATOMS,:)=0.0D0 
1392:               HESS(:, DEGFREEDOMS+1:3*NATOMS)=0.0D0 
1393:               HESS(1:DEGFREEDOMS, 1:DEGFREEDOMS)=XRIGIDHESS(1:DEGFREEDOMS, 1:DEGFREEDOMS) 
1394:            ENDIF 
1395:  
1396:        END IF 
1397:  
1398:        write(*,*) "Energy ", EREAL 
1399:        write(*,*) GRAD(:3*NATOMS) 
1400:  
1401:    END IF1371:    END IF
1402: !1372: !
1403: !  --------------- End of possible potentials - now add fields if required ------------------------------1373: !  --------------- End of possible potentials - now add fields if required ------------------------------
1404: !1374: !
1405: 1375: 
1406:    IF (RIGIDCONTOURT) THEN1376:    IF (RIGIDCONTOURT) THEN
1407: !       sf344> sanity check, since this keyword works only for a pair of rigid particles1377: !       sf344> sanity check, since this keyword works only for a pair of rigid particles
1408:       IF (NATOMS.NE.4) THEN1378:       IF (NATOMS.NE.4) THEN
1409:          WRITE(MYUNIT,'(A)') 'potential> ERROR - RIGIDCONTOUR works only for a PAIR of rigid bodies, exiting'1379:          WRITE(MYUNIT,'(A)') 'potential> ERROR - RIGIDCONTOUR works only for a PAIR of rigid bodies, exiting'
1410:          STOP1380:          STOP
1555: !            GRAD(J2-2)=GRAD(J2-2)-XG1525: !            GRAD(J2-2)=GRAD(J2-2)-XG
1556: !            GRAD(J2-1)=GRAD(J2-1)-YG1526: !            GRAD(J2-1)=GRAD(J2-1)-YG
1557: !            GRAD(J2)=  GRAD(J2)-ZG1527: !            GRAD(J2)=  GRAD(J2)-ZG
1558: !         END DO1528: !         END DO
1559: !      END IF1529: !      END IF
1560: 1530: 
1561:       IF (CSMT .AND. (.NOT.SYMMETRIZECSM)) THEN1531:       IF (CSMT .AND. (.NOT.SYMMETRIZECSM)) THEN
1562:          DUMMY2=0.0D01532:          DUMMY2=0.0D0
1563:          RMS=0.0D01533:          RMS=0.0D0
1564:       ELSE IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN1534:       ELSE IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
1565:          write(*,*) "Rigid body RMS" 
1566:          DUMMY2=SUM(GRAD(1:DEGFREEDOMS)**2)1535:          DUMMY2=SUM(GRAD(1:DEGFREEDOMS)**2)
1567:          RMS=MAX(SQRT(DUMMY2/DEGFREEDOMS), 1.0D-100)1536:          RMS=MAX(SQRT(DUMMY2/DEGFREEDOMS), 1.0D-100)
1568:          write(*,*) "RMS:", RMS 
1569:          IF ((RMS < 5.0D0 * BQMAX) .AND. AACONVERGENCET .AND. (.NOT. COMPRESSRIGIDT) .AND. (.NOT. COMPON)) THEN1537:          IF ((RMS < 5.0D0 * BQMAX) .AND. AACONVERGENCET .AND. (.NOT. COMPRESSRIGIDT) .AND. (.NOT. COMPON)) THEN
1570:             write(*,*) "Using AACONVERGENCE" 
1571:             CALL AACONVERGENCE (GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD, RMS)1538:             CALL AACONVERGENCE (GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD, RMS)
1572:             write(*,*) "RMS", RMS 
1573:          END IF1539:          END IF
1574:       ELSE IF (MLP3T.OR.MLPB3T.OR.MLQT) THEN1540:       ELSE IF (MLP3T.OR.MLPB3T.OR.MLQT) THEN
1575:          DUMMY2=SUM(GRAD(1:NATOMS)**2)1541:          DUMMY2=SUM(GRAD(1:NATOMS)**2)
1576:          RMS=MAX(DSQRT(DUMMY2/(NATOMS)), 1.0D-100)1542:          RMS=MAX(DSQRT(DUMMY2/(NATOMS)), 1.0D-100)
1577:       ELSE IF (.NOT.THOMSONT) THEN1543:       ELSE IF (.NOT.THOMSONT) THEN
1578:          DUMMY2=SUM(GRAD(1:3*NATOMS)**2)1544:          DUMMY2=SUM(GRAD(1:3*NATOMS)**2)
1579:          RMS=MAX(DSQRT(DUMMY2/(3*NATOMS)), 1.0D-100)1545:          RMS=MAX(DSQRT(DUMMY2/(3*NATOMS)), 1.0D-100)
1580:          write(*,*) "Using regular RMS", RMS 
1581:       ELSE1546:       ELSE
1582:          DUMMY2=SUM(GRAD(1:2*NATOMS)**2)1547:          DUMMY2=SUM(GRAD(1:2*NATOMS)**2)
1583:          RMS=MAX(DSQRT(DUMMY2/(2*NATOMS)), 1.0D-100)1548:          RMS=MAX(DSQRT(DUMMY2/(2*NATOMS)), 1.0D-100)
1584:       END IF1549:       END IF
1585:       IF(DEBUG.AND.(RMS.NE.RMS)) THEN1550:       IF(DEBUG.AND.(RMS.NE.RMS)) THEN
1586:          WRITE(MYUNIT,'(A)' ) 'potential> WARNING - RMS force is NaN - if using AMBER igb=1, can be due to negative Born radii'1551:          WRITE(MYUNIT,'(A)' ) 'potential> WARNING - RMS force is NaN - if using AMBER igb=1, can be due to negative Born radii'
1587:       END IF1552:       END IF
1588: 1553: 
1589: !   IF ((CSMRMS < GUIDECUT).AND.(CSMDOGUIDET)) THEN1554: !   IF ((CSMRMS < GUIDECUT).AND.(CSMDOGUIDET)) THEN
1590:       IF (CSMDOGUIDET) THEN1555:       IF (CSMDOGUIDET) THEN


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0