hdiff output

r29298/potential.f 2015-11-17 23:32:51.596392044 +0000 r29297/potential.f 2015-11-17 23:32:52.192400037 +0000
3012:                   XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)3012:                   XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)
3013:                   CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)3013:                   CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)
3014:                ENDIF3014:                ENDIF
3015:                ! hk286 - numerical second derivatives3015:                ! hk286 - numerical second derivatives
3016:                ! CALL NABSECDER(COORDS,STEST)3016:                ! CALL NABSECDER(COORDS,STEST)
3017: 3017: 
3018:                IF (ALLOCATED(HESS)) DEALLOCATE(HESS)3018:                IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
3019:                IF (.NOT.ALLOCATED(TEMPHESS)) ALLOCATE(TEMPHESS(9*NATOMS*NATOMS))3019:                IF (.NOT.ALLOCATED(TEMPHESS)) ALLOCATE(TEMPHESS(9*NATOMS*NATOMS))
3020:                TEMPHESS(:) = 0.0D03020:                TEMPHESS(:) = 0.0D0
3021:                CALL MME2WRAPPER(COORDS,ENERGY,GRADATOMS,TEMPHESS,ATMASS,GRAD1)3021:                CALL MME2WRAPPER(COORDS,ENERGY,GRADATOMS,TEMPHESS,ATMASS,GRAD1)
3022:                IF (SQRT(ABS(ENERGY)) > 1.0D40) THEN 
3023:                   WRITE(*, *) "Linear dihedral detected in NAB routines (sff2.c)." 
3024:                END IF 
3025:                VNEW(1:3*NATOMS) = GRADATOMS(:)3022:                VNEW(1:3*NATOMS) = GRADATOMS(:)
3026:                ALLOCATE(HESS(3*NATOMS,3*NATOMS))3023:                ALLOCATE(HESS(3*NATOMS,3*NATOMS))
3027:                k=13024:                k=1
3028:                DO i=1,3*NATOMS3025:                DO i=1,3*NATOMS
3029:                   DO j=1,3*NATOMS3026:                   DO j=1,3*NATOMS
3030:                      HESS(i,j) = TEMPHESS(k)3027:                      HESS(i,j) = TEMPHESS(k)
3031:                      k=k+13028:                      k=k+1
3032:                   END DO3029:                   END DO
3033:                END DO3030:                END DO
3034:                DEALLOCATE(TEMPHESS)3031:                DEALLOCATE(TEMPHESS)


r29298/potential.f90 2015-11-17 23:32:51.212386894 +0000 r29297/potential.f90 2015-11-17 23:32:51.788394619 +0000
434:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)               434:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)               
435:          END IF435:          END IF
436:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)436:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
437:          ALLOCATE(HESS(3*NATOMS, 3*NATOMS))437:          ALLOCATE(HESS(3*NATOMS, 3*NATOMS))
438:          IF (CUDAT) THEN438:          IF (CUDAT) THEN
439:             CALL CUDA_NUMERICAL_HESS(NATOMS, X, HESS, DELTA=1.0D-6)439:             CALL CUDA_NUMERICAL_HESS(NATOMS, X, HESS, DELTA=1.0D-6)
440:          ELSE 440:          ELSE 
441:             IF (.NOT. ALLOCATED(TEMPHESS)) ALLOCATE(TEMPHESS(9*NATOMS*NATOMS))441:             IF (.NOT. ALLOCATED(TEMPHESS)) ALLOCATE(TEMPHESS(9*NATOMS*NATOMS))
442:             TEMPHESS(:)=0.0D0442:             TEMPHESS(:)=0.0D0
443:             CALL MME2WRAPPER(X, ENERGY, VNEW, TEMPHESS, ATMASS1, GRAD1)443:             CALL MME2WRAPPER(X, ENERGY, VNEW, TEMPHESS, ATMASS1, GRAD1)
444:             IF (SQRT(ABS(ENERGY)) > 1.0D40) THEN 
445:                EREAL = ENERGY 
446:                WRITE(MYUNIT, *) "Linear dihedral detected in NAB routines (sff2.c)." 
447:             END IF 
448:             K=1444:             K=1
449:             DO I=1, 3*NATOMS445:             DO I=1, 3*NATOMS
450:                DO J=1, 3*NATOMS446:                DO J=1, 3*NATOMS
451:                   HESS(I, J)=TEMPHESS(K)447:                   HESS(I, J)=TEMPHESS(K)
452:                   K=K + 1448:                   K=K + 1
453:                END DO449:                END DO
454:             END DO450:             END DO
455:             DEALLOCATE(TEMPHESS)451:             DEALLOCATE(TEMPHESS)
456:          END IF452:          END IF
457:          IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN453:          IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
553:          VNEW=0.0D0549:          VNEW=0.0D0
554: ! khs26> Copied analytical second derivatives from OPTIM550: ! khs26> Copied analytical second derivatives from OPTIM
555:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN551:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
556:             XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)552:             XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
557:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)               553:             CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, X, XRIGIDCOORDS)               
558:          END IF554:          END IF
559:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)555:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
560:          IF (.NOT. ALLOCATED(TEMPHESS)) ALLOCATE(TEMPHESS(9*NATOMS*NATOMS))556:          IF (.NOT. ALLOCATED(TEMPHESS)) ALLOCATE(TEMPHESS(9*NATOMS*NATOMS))
561:          TEMPHESS(:)=0.0D0557:          TEMPHESS(:)=0.0D0
562:          CALL MME2WRAPPER(X, ENERGY, VNEW, TEMPHESS, ATMASS1, GRAD1)558:          CALL MME2WRAPPER(X, ENERGY, VNEW, TEMPHESS, ATMASS1, GRAD1)
563:          IF (SQRT(ABS(ENERGY)) > 1.0D40) THEN 
564:             EREAL = ENERGY 
565:             WRITE(MYUNIT, *) "Linear dihedral detected in NAB routines (sff2.c)." 
566:          END IF 
567:          ALLOCATE(HESS(3*NATOMS, 3*NATOMS))559:          ALLOCATE(HESS(3*NATOMS, 3*NATOMS))
568:          K=1560:          K=1
569:          DO I=1, 3*NATOMS561:          DO I=1, 3*NATOMS
570:             DO J=1, 3*NATOMS562:             DO J=1, 3*NATOMS
571:                HESS(I, J)=TEMPHESS(K)563:                HESS(I, J)=TEMPHESS(K)
572:                K=K + 1564:                K=K + 1
573:             END DO565:             END DO
574:          END DO566:          END DO
575:          DEALLOCATE(TEMPHESS)567:          DEALLOCATE(TEMPHESS)
576:          IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN568:          IF ( RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN


r29298/sff2.c 2015-11-17 23:32:51.404389468 +0000 r29297/sff2.c 2015-11-17 23:32:52.000397462 +0000
1246:     gz = ykj * xkl - xkj * ykl;1246:     gz = ykj * xkl - xkj * ykl;
1247: 1247: 
1248:     bi = dx * dx + dy * dy + dz * dz;1248:     bi = dx * dx + dy * dy + dz * dz;
1249:     bk = gx * gx + gy * gy + gz * gz;1249:     bk = gx * gx + gy * gy + gz * gz;
1250:     ct = dx * gx + dy * gy + dz * gz;1250:     ct = dx * gx + dy * gy + dz * gz;
1251: 1251: 
1252:       s = xkj*(dz*gy-dy*gz)+ykj*(dx*gz-dz*gx)+ zkj*(dy*gx-dx*gy);1252:       s = xkj*(dz*gy-dy*gz)+ykj*(dx*gz-dz*gx)+ zkj*(dy*gx-dx*gy);
1253: 1253: 
1254:     /*    ----- branch if linear dihedral -----   */1254:     /*    ----- branch if linear dihedral -----   */
1255: 1255: 
1256:     if (bi < 0.01 || bk < 0.01) { 
1257:        /* Set the energy to be large, such that it isn't used. */ 
1258:        e_tors = 1.0e100; 
1259:        /* Set the hessian to be the identity matrix. */ 
1260:        /* N.B. n = 3*natoms */ 
1261:             for (i = 0; i < n; i++) { 
1262:           for (j = 0; j < n; j++) { 
1263:              g[j + n * i] = 0.0; 
1264:           } 
1265:        } 
1266:        for (i = 0; i < n; i++) { 
1267:           g[i + n * i] = 1.0; 
1268:        } 
1269:        return e_tors; 
1270:     } 
1271:     /* We used to use these assertions to handle linear dihedrals. 
1272:        This causes GMIN to terminate if a linear dihedral was encountered. 
1273:     assert(bi > 0.01);1256:     assert(bi > 0.01);
1274:     assert(bk > 0.01);1257:     assert(bk > 0.01);
1275:     */ 
1276: 1258: 
1277:     boi2 = 1. / bi;1259:     boi2 = 1. / bi;
1278:     boj2 = 1. / bk;1260:     boj2 = 1. / bk;
1279:     bi = sqrt(bi);1261:     bi = sqrt(bi);
1280:     bk = sqrt(bk);1262:     bk = sqrt(bk);
1281:     z1 = 1. / bi;1263:     z1 = 1. / bi;
1282:     z2 = 1. / bk;1264:     z2 = 1. / bk;
1283:     bioj = bi * z2;1265:     bioj = bi * z2;
1284:     bjoi = bk * z1;1266:     bjoi = bk * z1;
1285:     ct = ct * z1 * z2;1267:     ct = ct * z1 * z2;
1575: #endif1557: #endif
1576: 1558: 
1577:     /*1559:     /*
1578:      *     ----- now set up the ddr array = second der. of cosphi w/respect1560:      *     ----- now set up the ddr array = second der. of cosphi w/respect
1579:      *           to cartesians; first we take the first term of last formula1561:      *           to cartesians; first we take the first term of last formula
1580:      *           on p. 113 of cff book -----1562:      *           on p. 113 of cff book -----
1581:      */1563:      */
1582: 1564: 
1583:     for (i = 1; i <= 12; i++) {1565:     for (i = 1; i <= 12; i++) {
1584:       for (j = i; j <= 12; j++) {1566:       for (j = i; j <= 12; j++) {
1585:              ddr[i][j] = 0.0;1567:         ddr[i][j] = 0.0;
1586:              for (k = 1; k <= 6; k++) {1568:         for (k = 1; k <= 6; k++) {
1587:                for (l = 1; l <= 6; l++) {1569:           for (l = 1; l <= 6; l++) {
1588:                  ddr[i][j] = ddr[i][j] + ddc[k][l]*dtx[k][i]*dtx[l][j];1570:             ddr[i][j] = 
1589:                }1571:               ddr[i][j] + ddc[k][l]*dtx[k][i]*dtx[l][j];
1590:              }1572:           }
 1573:         }
1591:       }1574:       }
1592:     }1575:     }
1593: 1576: 
1594:     /*     ----- now do the second term of this equation -----  */1577:     /*     ----- now do the second term of this equation -----  */
1595: 1578: 
1596:     ddr[2][9] += dc[1];1579:     ddr[2][9] += dc[1];
1597:     ddr[3][5] += dc[1];1580:     ddr[3][5] += dc[1];
1598:     ddr[6][8] += dc[1];1581:     ddr[6][8] += dc[1];
1599:     ddr[2][6] -= dc[1];1582:     ddr[2][6] -= dc[1];
1600:     ddr[3][8] -= dc[1];1583:     ddr[3][8] -= dc[1];
1633:     /*     ----- Now form the second derivative matrix -----  */1616:     /*     ----- Now form the second derivative matrix -----  */
1634: 1617: 
1635:     for (i = 1; i <= 12; i++) {1618:     for (i = 1; i <= 12; i++) {
1636:       ist = at4;1619:       ist = at4;
1637:       if (i <= 9) ist = at3;1620:       if (i <= 9) ist = at3;
1638:       if (i <= 6) ist = at2;1621:       if (i <= 6) ist = at2;
1639:       if (i <= 3) ist = at1;1622:       if (i <= 3) ist = at1;
1640:       iof = i % 3; if (iof == 0) iof = 3;1623:       iof = i % 3; if (iof == 0) iof = 3;
1641:       inew = ist + iof - 1;1624:       inew = ist + iof - 1;
1642:       for (j = i; j <= 12; j++) {1625:       for (j = i; j <= 12; j++) {
1643:              jst = at4;1626:         jst = at4;
1644:         if (j <= 9) jst = at3;1627:         if (j <= 9) jst = at3;
1645:         if (j <= 6) jst = at2;1628:         if (j <= 6) jst = at2;
1646:         if (j <= 3) jst = at1;1629:         if (j <= 3) jst = at1;
1647:         jof = j % 3; if (jof == 0) jof = 3;1630:         jof = j % 3; if (jof == 0) jof = 3;
1648:         jnew = jst + jof - 1;1631:         jnew = jst + jof - 1;
 1632: 
1649: #ifdef SCALAPACK1633: #ifdef SCALAPACK
1650:         ptr = ptr2d(g, descG_PxQ, jnew, inew);1634: 
1651:         if (ptr != NULL) *ptr += ddf*dr[i]*dr[j] + df*ddr[i][j];1635:         ptr = ptr2d(g, descG_PxQ, jnew, inew);
1652:         if (inew != jnew) {1636:         if (ptr != NULL) *ptr += ddf*dr[i]*dr[j] + df*ddr[i][j];
1653:           ptr = ptr2d(g, descG_PxQ, inew, jnew);1637:         if (inew != jnew) {
1654:           if (ptr != NULL) *ptr += ddf*dr[i]*dr[j] + df*ddr[i][j];1638:           ptr = ptr2d(g, descG_PxQ, inew, jnew);
1655:              }1639:           if (ptr != NULL) *ptr += ddf*dr[i]*dr[j] + df*ddr[i][j];
 1640:         }
1656: #else1641: #else
1657:              g[inew + n*jnew] += ddf*dr[i]*dr[j] + df*ddr[i][j];1642:         g[inew + n*jnew] += ddf*dr[i]*dr[j] + df*ddr[i][j];
1658:              if (inew != jnew)1643:         if (inew != jnew)
1659:                g[jnew + n*inew] += ddf*dr[i]*dr[j] + df*ddr[i][j];1644:           g[jnew + n*inew] += ddf*dr[i]*dr[j] + df*ddr[i][j];
1660: #endif1645: #endif
1661:       }1646:       }
1662:     }1647:     }
1663: 1648: 
1664: 1649: 
1665: #ifdef PRINT_EPHI1650: #ifdef PRINT_EPHI
1666:     printf("%4d %4d %4d %4d %4d %9.4f\n", i + 1, at1 / 3, at2 / 3,1651:     printf("%4d %4d %4d %4d %4d %9.4f\n", i + 1, at1 / 3, at2 / 3,
1667:            at3 / 3, at4 / 3, e);1652:            at3 / 3, at4 / 3, e);
1668: #endif1653: #endif
1669: 1654: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0