hdiff output

r14089/amber9dummy.f90 2017-01-21 10:32:30.551188235 +0000 r14088/amber9dummy.f90 2017-01-21 10:32:35.031188235 +0000
112: subroutine GETICCOORDS()112: subroutine GETICCOORDS()
113: 113: 
114: end subroutine GETICCOORDS114: end subroutine GETICCOORDS
115: 115: 
116: subroutine CHECK_SIDECHAIN(i3,j3,k3,l3,IICD,IS_SIDECHAIN)116: subroutine CHECK_SIDECHAIN(i3,j3,k3,l3,IICD,IS_SIDECHAIN)
117: integer :: i3,j3,k3,l3,IICD117: integer :: i3,j3,k3,l3,IICD
118: logical :: IS_SIDECHAIN(100)118: logical :: IS_SIDECHAIN(100)
119: 119: 
120: end subroutine CHECK_SIDECHAIN120: end subroutine CHECK_SIDECHAIN
121: 121: 
122: subroutine CHGETICVAL(Q,BOND1, BOND2, THET1, THET2, PHI, DIHED_ONLY)122: subroutine CHGETICVAL(Q,BOND1, BOND2, THET1, THET2, PHI)
123: double precision :: Q(3*1000), BOND1(1000), BOND2(1000), THET1(1000)123: double precision :: Q(3*1000), BOND1(1000), BOND2(1000), THET1(1000)
124: double precision :: THET2(1000), PHI(1000)124: double precision :: THET2(1000), PHI(1000)
125: logical :: DIHED_ONLY 
126: end subroutine CHGETICVAL125: end subroutine CHGETICVAL
127: 126: 
128: double precision function PHI_INTERP(ST_PHI,FI_PHI,SFRAC) 
129: double precision:: ST_PHI,FI_PHI,PHI,SFRAC 
130: end function PHI_INTERP 
131:  
132: subroutine  TAKESTEPAMDIHED(COORDS, CSTART, CFINISH,SFRAC)127: subroutine  TAKESTEPAMDIHED(COORDS, CSTART, CFINISH,SFRAC)
133: double precision :: COORDS(3*1000), CSTART(3*1000), CFINISH(3*1000), SFRAC128: double precision :: COORDS(3*1000), CSTART(3*1000), CFINISH(3*1000), SFRAC
134: 129: 
135: end subroutine TAKESTEPAMDIHED130: end subroutine TAKESTEPAMDIHED


r14089/amberinterface.f 2017-01-21 10:32:27.827188235 +0000 r14088/amberinterface.f 2017-01-21 10:32:32.635188235 +0000
1354:       double precision :: rij, xji, ant , rij2,ant21354:       double precision :: rij, xji, ant , rij2,ant2
1355:       double precision :: vec_ij(3), vec_kj(3)1355:       double precision :: vec_ij(3), vec_kj(3)
1356:       integer :: i,j, i3, j3, k3, l3, count1356:       integer :: i,j, i3, j3, k3, l3, count
1357: !      INTEGER :: IC_COORDS(nres*95+5*nres,4)1357: !      INTEGER :: IC_COORDS(nres*95+5*nres,4)
1358: !      CHARACTER(LEN=3) :: IC_IMPROP(nres*95+5*nres)1358: !      CHARACTER(LEN=3) :: IC_IMPROP(nres*95+5*nres)
1359: !      INTEGER :: lenic1359: !      INTEGER :: lenic
1360:       double precision pt9991360:       double precision pt999
1361:       data pt999 /0.9990d0/1361:       data pt999 /0.9990d0/
1362: !ask whether you have to consider the qpsander keyword 1362: !ask whether you have to consider the qpsander keyword 
1363: 1363: 
 1364:       PRINT*, "in geopt"
1364:       count = 01365:       count = 0
1365:       DO i = 0, nphih-11366:       DO i = 0, nphih-1
1366:         i3 = ix(i40+i)/3 + 11367:         i3 = ix(i40+i)/3 + 1
1367:         j3 = ix(i42+i)/3 + 11368:         j3 = ix(i42+i)/3 + 1
1368:         k3 = ix(i44+i)1369:         k3 = ix(i44+i)
1369: ! NOTE: there is a BUG for rows with negative JUST on 4th term1370: ! NOTE: there is a BUG for rows with negative JUST on 4th term
1370: !                         if they exist1371: !                         if they exist
1371:         l3 = ix(i46+i)1372:         l3 = ix(i46+i)
1372:         IF (k3*l3>0) THEN1373:         IF (k3*l3>0) THEN
1373:            count = count + 11374:            count = count + 1
1464:       ENDDO1465:       ENDDO
1465:       END SUBROUTINE FINDALLCOORDS1466:       END SUBROUTINE FINDALLCOORDS
1466: 1467: 
1467: 1468: 
1468: ! ***********************************************************************1469: ! ***********************************************************************
1469: 1470: 
1470: 1471: 
1471:       SUBROUTINE TAKESTEPAMDIHED(Q, CSTART, CFINISH, SFRAC)1472:       SUBROUTINE TAKESTEPAMDIHED(Q, CSTART, CFINISH, SFRAC)
1472:       use modamber91473:       use modamber9
1473:       use commons, only : natoms1474:       use commons, only : natoms
1474:       IMPLICIT NONE 
1475: ! similar to ICINTERP in twist.src1475: ! similar to ICINTERP in twist.src
1476: ! BOND1, THET1, PHI, THET2, BOND2 - IC arrays for interpolated structure1476: ! BOND1, THET1, PHI, THET2, BOND2 - IC arrays for interpolated structure
1477: ! GETICCOORDS generates the atoms defining the internal coordinates1477: ! GETICCOORDS generates the atoms defining the internal coordinates
1478: ! CHGETICVAL fills the internal coordinate table (similar to FILLICTABLE in Charmm)1478: ! CHGETICVAL fills the internal coordinate table (similar to FILLICTABLE in Charmm)
1479:       DOUBLE PRECISION, INTENT(IN) :: CSTART(3*natoms),CFINISH(3*natoms)1479:       DOUBLE PRECISION, INTENT(IN) :: CSTART(3*natoms),CFINISH(3*natoms)
1480:       DOUBLE PRECISION, INTENT(INOUT) :: Q(3*natoms)1480:       DOUBLE PRECISION, INTENT(INOUT) :: Q(3*natoms)
1481:       DOUBLE PRECISION:: TEMP_AR(9)1481:       DOUBLE PRECISION:: TEMP_AR(9)
1482:       DOUBLE PRECISION :: PHIINTERP 
1483:       DOUBLE PRECISION :: SFRAC1482:       DOUBLE PRECISION :: SFRAC
1484:       double precision :: ST_PHI(nphia+nphih), FI_PHI(nphia+nphih), &1483:       double precision :: ST_PHI(nphia+nphih), FI_PHI(nphia+nphih), &
1485:      &      ST_THET1(ntheth+ntheta+ngper),ST_THET2(ntheth+ntheta+ngper), &1484:      &      ST_THET1(ntheth+ntheta+ngper),ST_THET2(ntheth+ntheta+ngper), &
1486:      &      FI_THET1(ntheth+ntheta+ngper),FI_THET2(ntheth+ntheta+ngper), &1485:      &      FI_THET1(ntheth+ntheta+ngper),FI_THET2(ntheth+ntheta+ngper), &
1487:      &      PHI(nphia+nphih), THET1(ntheth+ntheta+ngper), &1486:      &      PHI(nphia+nphih), THET1(ntheth+ntheta+ngper), &
1488:      &      THET2(ntheth+ntheta+ngper)1487:      &      THET2(ntheth+ntheta+ngper)
1489:       double precision :: ST_BOND1(nbonh+nbona+nbper), &1488:       double precision :: ST_BOND1(nbonh+nbona+nbper), &
1490:      &      ST_BOND2(nbonh+nbona+nbper), FI_BOND1(nbonh+nbona+nbper),&1489:      &      ST_BOND2(nbonh+nbona+nbper), FI_BOND1(nbonh+nbona+nbper),&
1491:      &      FI_BOND2(nbonh+nbona+nbper), &1490:      &      FI_BOND2(nbonh+nbona+nbper), &
1492:      &      BOND1(nbonh+nbona+nbper),BOND2(nbonh+nbona+nbper) 1491:      &      BOND1(nbonh+nbona+nbper),BOND2(nbonh+nbona+nbper) 
1493:       INTEGER :: i, nterm_res, lenic2, IICD,I1,J11492: !      INTEGER :: IC_COORDS(nres*95+5*nres,4)
 1493: !      CHARACTER(LEN=3) :: IC_IMPROP(nres*95+5*nres)
 1494:       INTEGER :: i, nterm_res, lenic2
1494:       INTEGER :: atomOffset, nat_res,k, at1, at2,at31495:       INTEGER :: atomOffset, nat_res,k, at1, at2,at3
1495:       LOGICAL :: found, found1, found2, found3, BACKINTERP1496:       LOGICAL :: found, found1, found2, found3
1496:       DOUBLE PRECISION :: ST_NCOOR(3*natoms), FI_NCOOR(3*natoms)1497:       DOUBLE PRECISION :: ST_NCOOR(3*natoms), FI_NCOOR(3*natoms)
1497:       double precision ::pt9991498:       double precision pt999
1498:       DOUBLE PRECISION :: PI, RADDEG, DEGRAD 
1499:       data pt999 /0.9990d0/1499:       data pt999 /0.9990d0/
1500:       PARAMETER(PI=3.141592653589793D0) 
1501:       PARAMETER (RADDEG=180.0D0/PI) 
1502:       PARAMETER (DEGRAD=PI/180.0D0) 
1503: !ask whether you have to consider the qpsander keyword 1500: !ask whether you have to consider the qpsander keyword 
1504: 1501: 
1505:       nterm_res =ix(i02+1)- ix(i02)1502:       nterm_res =ix(i02+1)- ix(i02)
1506: 1503: 
1507: 1504: 
1508:       PRINT*, "IC"1505:       PRINT*, "IC"
1509:       IF (.NOT. ALLOCATED(IC_COORDS)) THEN1506:       IF (.NOT. ALLOCATED(IC_COORDS)) THEN
1510:          CALL GETICCOORDS()1507:          CALL GETICCOORDS()
1511:          PRINT*, "lenic", lenic1508:          PRINT*, "lenic", lenic
1512:       ENDIF1509:       ENDIF
1513: 1510: 
 1511:       lenic2 = lenic
 1512:       
1514:       IF (.NOT. ALLOCATED(IS_SIDECHAIN)) THEN1513:       IF (.NOT. ALLOCATED(IS_SIDECHAIN)) THEN
1515:           ALLOCATE(IS_SIDECHAIN(lenic))1514:           ALLOCATE(IS_SIDECHAIN(lenic))
1516:           DO i = 1, lenic1515:           DO i = 1, lenic
1517:               CALL CHECK_SIDECHAIN(IC_COORDS(i,1), IC_COORDS(i,2),      &1516:               CALL CHECK_SIDECHAIN(IC_COORDS(i,1), IC_COORDS(i,2),&
1518:      &        IC_COORDS(i,3),IC_COORDS(i,4), i, IS_SIDECHAIN)1517:      &        IC_COORDS(i,3),IC_COORDS(i,4), i, IS_SIDECHAIN)
1519:           ENDDO1518:           ENDDO
1520:       ENDIF1519:       ENDIF
1521: 1520:  
1522:       lenic2 = lenic 
1523:        
1524:       CALL CHGETICVAL(CSTART, ST_BOND1, &1521:       CALL CHGETICVAL(CSTART, ST_BOND1, &
1525:      &     ST_BOND2,ST_THET1, ST_THET2, ST_PHI,.FALSE.)1522:      &     ST_BOND2,ST_THET1, ST_THET2, ST_PHI)
1526:       CALL CHGETICVAL(CFINISH, FI_BOND1, &1523:       CALL CHGETICVAL(CFINISH, FI_BOND1, &
1527:      &    FI_BOND2, FI_THET1, FI_THET2, FI_PHI,.FALSE.)1524:      &    FI_BOND2, FI_THET1, FI_THET2, FI_PHI)
1528:       1525:       
1529: !      DO i = 1, lenic1526: ! the following 14lines are just a test 
1530: !         PRINT '(i4,4a5,2f7.4,2f9.3,2f10.3)', i,&1527: !      ST_NCOOR =9999.0000000
1531: !     &    ih(m04+IC_COORDS(i,1)-1),ih(m04+IC_COORDS(i,2)-1),&1528: !first three coords - N,CA,C
1532: !     &    ih(m04+IC_COORDS(i,3)-1),ih(m04+IC_COORDS(i,4)-1),&1529: !      CALL GETSEEDATM("N","CA","C", at1, at2,at3)
1533: !     &    ST_BOND1(i),ST_BOND2(i),ST_THET1(i),ST_THET2(i),ST_PHI(i), FI_PHI(i)1530: !      PRINT*, "ats", at1, at2, at3
1534: !      ENDDO 
1535:  
1536:       IF (AMBPERTT) THEN 
1537:          CALL GET_TWISTABLE(ST_PHI,FI_PHI) 
1538:       ENDIF 
1539:  
1540: !      CALL CH_SEED(at1,at2,at3,ST_NCOOR,ST_BOND1, ST_BOND2,1531: !      CALL CH_SEED(at1,at2,at3,ST_NCOOR,ST_BOND1, ST_BOND2,
1541: !     &     ST_THET1,ST_THET2, ST_PHI)1532: !     &     ST_THET1,ST_THET2, ST_PHI)
 1533: !      CALL CHARMMBILDC(1, lenic2, ST_NCOOR, ST_BOND1,ST_BOND2,ST_THET1,
 1534: !     &     ST_THET2, ST_PHI)
 1535: !      DO i =1, NATOMS
 1536: !         PRINT '(a4,4f11.5)',ih(m04+i-1),ST_NCOOR(3*(i-1)+1),
 1537: !     &          ST_NCOOR(3*(i-1)+2),ST_NCOOR(3*(i-1)+3)
 1538: !      ENDDO
 1539: ! end test
1542: 1540: 
 1541:       CALL CHGETICVAL(Q, BOND1,   BOND2, THET1, THET2, PHI)
 1542: !      DO i=1, lenic
 1543: !         PRINT '(i4,l,4a5,f7.4,3f10.3,f7.4)', i, IS_SIDECHAIN(i),
 1544: !     &    ih(m04+IC_COORDS(i,1)-1),
 1545: !     &    ih(m04+IC_COORDS(i,2)-1), ih(m04+IC_COORDS(i,3)-1),
 1546: !     &    ih(m04+IC_COORDS(i,4)-1), BOND1(i), THET1(i), PHI(i),
 1547: !     &    THET2(i), BOND2(i)
 1548: !      ENDDO
1543: 1549: 
1544:       CALL CHGETICVAL(Q, BOND1, BOND2, THET1, THET2, PHI,.FALSE.)1550:       DO i =1,lenic
1545: 1551:          IF (IS_SIDECHAIN(i) .EQV. .TRUE.) THEN
1546:       I1 = 11552:             BOND1(i) = SFRAC*ST_BOND1(i)+(1.D0-SFRAC)*FI_BOND1(i)
1547:       J1 = 11553:             BOND2(i) = SFRAC*ST_BOND2(i)+(1.D0-SFRAC)*FI_BOND2(i)
1548:       DO IICD =1,lenic1554:             THET1(i) = SFRAC*ST_THET1(i)+(1.D0-SFRAC)*FI_THET1(i)
1549:           BACKINTERP = .FALSE.1555:             THET2(i) = SFRAC*ST_THET2(i)+(1.D0-SFRAC)*FI_THET2(i)
1550:           IF (AMBIT) THEN !backbone interpolation in internals too1556:             IF (ST_PHI(i) .GE. 0.D0 .AND. FI_PHI(i) .GE. 0.D0) THEN 
1551: !            Do not interpolate bond lengths and angles;1557:                PHI(i) = SFRAC*ST_PHI(i) + (1.D0-SFRAC)*FI_PHI(i)
1552: !            instead take the value from the closer endpoint.1558:             ELSEIF (ST_PHI(i) .LT. 0.D0 .AND. FI_PHI(i) .LT. 0.D0) THEN
1553:              IF (SFRAC.GE.0.5D0) THEN1559:                PHI(i) = SFRAC*ST_PHI(i) + (1.D0-SFRAC)*FI_PHI(i)
1554:                   BOND1(IICD)=ST_BOND1(IICD)1560:             ELSEIF (ST_PHI(i) .GE. 0.D0 .AND. FI_PHI(i) .LT. 0.D0) THEN
1555:                   BOND2(IICD)=ST_BOND2(IICD)1561:                IF (DABS(ST_PHI(i)-FI_PHI(i)).LE.180.D0) THEN
1556:                   THET1(IICD)=ST_THET1(IICD)1562:                   PHI(i) = SFRAC*ST_PHI(i) + (1.D0-SFRAC)*FI_PHI(i)
1557:                   THET2(IICD)=ST_THET2(IICD)1563:                ELSE
1558:              ELSE1564:                   IF (ST_PHI(i).GT.DABS(FI_PHI(i))) THEN
1559:                   BOND1(IICD)=FI_BOND1(IICD)1565:                      ST_PHI(i)=ST_PHI(i)-360.D0
1560:                   BOND2(IICD)=FI_BOND2(IICD)1566:                   ELSE
1561:                   THET1(IICD)=FI_THET1(IICD)1567:                      ST_PHI(i)=FI_PHI(i)+360.D0
1562:                   THET2(IICD)=FI_THET2(IICD)1568:                   ENDIF
1563:              ENDIF1569:                   PHI(i) = SFRAC*ST_PHI(i) + (1.D0-SFRAC)*FI_PHI(i)
1564:              IF (IICD .EQ. PHIPSI(J1)) THEN !interpolate Ramachandran angles1570:                ENDIF
1565:                  J1 = J1 +11571:             ELSEIF (ST_PHI(i) .LT. 0.D0 .AND. FI_PHI(i) .GE. 0.D0) THEN
1566:                  PHI(IICD) = PHIINTERP(ST_PHI(IICD),FI_PHI(IICD),SFRAC)1572:                IF (DABS(ST_PHI(i)-FI_PHI(i)).LE.180.D0) THEN
1567:              ELSEIF (IICD .EQ. TW_SIDECHAIN(I1)) THEN1573:                   PHI(i) = SFRAC*ST_PHI(i) + (1.D0-SFRAC)*FI_PHI(i)
1568:                  I1 = I1 +11574:                ELSE
1569:                  PHI(IICD) = PHIINTERP(ST_PHI(IICD), FI_PHI(IICD), SFRAC)1575:                   IF (DABS(ST_PHI(i)).GT.FI_PHI(i)) THEN
1570:              ELSE 1576:                      ST_PHI(i)=ST_PHI(i)+360.D0
1571:                  IF (SFRAC.GE.0.5D0) THEN1577:                   ELSE
1572:                     PHI(IICD)=ST_PHI(IICD)1578:                      FI_PHI(i)=FI_PHI(i)-360.D0
1573:                  ELSE1579:                   ENDIF
1574:                     PHI(IICD)=FI_PHI(IICD)1580:                   PHI(i) = SFRAC*ST_PHI(i) + (1.D0-SFRAC)*FI_PHI(i)
1575:                  ENDIF1581:                ENDIF
1576:              ENDIF    1582:             IF (PHI(i).GT.180.D0) PHI(i)=PHI(i)-360.D0
1577:           ELSE1583:             IF (PHI(i).LE.-180.D0) PHI(i)=PHI(i)+360.D0
1578:           1584:             ENDIF
1579: !            Twistable dihedrals are interpolated,1585:          ENDIF
1580: !            for nontwistable dihedrals and impropers 
1581: !            the value from the closer endpoint is used. 
1582:              IF (IS_SIDECHAIN(IICD) .EQV. .TRUE.) THEN !all sidechains 
1583:                  IF (SFRAC.GE.0.5D0) THEN 
1584:                       BOND1(IICD)=ST_BOND1(IICD) 
1585:                       BOND2(IICD)=ST_BOND2(IICD) 
1586:                       THET1(IICD)=ST_THET1(IICD) 
1587:                       THET2(IICD)=ST_THET2(IICD) 
1588:                  ELSE 
1589:                       BOND1(IICD)=FI_BOND1(IICD) 
1590:                       BOND2(IICD)=FI_BOND2(IICD) 
1591:                       THET1(IICD)=FI_THET1(IICD) 
1592:                       THET2(IICD)=FI_THET2(IICD) 
1593:                  ENDIF 
1594:              !take internal coord value for all sidechains to make them look decent 
1595:                  IF (IICD .EQ. TW_SIDECHAIN(I1)) THEN 
1596:                     I1 = I1 + 1 
1597:                     PHI(IICD) = PHIINTERP(ST_PHI(IICD), FI_PHI(IICD), SFRAC) 
1598:                  ELSE  
1599:                     IF (SFRAC.GE.0.5D0) THEN 
1600:                        PHI(IICD)=ST_PHI(IICD) 
1601:                     ELSE 
1602:                        PHI(IICD)=FI_PHI(IICD) 
1603:                     ENDIF 
1604:                  ENDIF 
1605:               ENDIF  
1606:           ENDIF 
1607:       ENDDO1586:       ENDDO
1608: 1587: 
1609: !      PRINT*, "interpolated"1588: !      PRINT*, "interpolated"
1610: !      DO i=1, lenic1589: !      DO i=1, lenic
1611: !         PRINT '(i4,l,4a5,f7.4,3f10.3,f7.4)', i, IS_SIDECHAIN(i),&1590: !         PRINT '(i4,l,4a5,f7.4,3f10.3,f7.4)', i, IS_SIDECHAIN(i),
1612: !     &    ih(m04+IC_COORDS(i,1)-1) ,                             &1591: !     &    ih(m04+IC_COORDS(i,1)-1),
1613: !     &    ih(m04+IC_COORDS(i,2)-1), ih(m04+IC_COORDS(i,3)-1),    &1592: !     &    ih(m04+IC_COORDS(i,2)-1), ih(m04+IC_COORDS(i,3)-1),
1614: !     &    ih(m04+IC_COORDS(i,4)-1), BOND1(i), THET1(i), PHI(i),  &1593: !     &    ih(m04+IC_COORDS(i,4)-1), BOND1(i), THET1(i), PHI(i),
1615: !     &    THET2(i), BOND2(i)1594: !     &    THET2(i), BOND2(i)
1616: !      ENDDO1595: !      ENDDO
1617: 1596: 
1618:       CALL GETSEEDATM("N","CA", "C", at1,at2,at3,0)1597:       CALL GETSEEDATM("N","CA", "C", at1,at2,at3,0)
1619:       TEMP_AR(1)= Q(3*(at1-1)+1);TEMP_AR(2)= Q(3*(at1-1)+2)1598:       TEMP_AR(1)= Q(3*(at1-1)+1);TEMP_AR(2)= Q(3*(at1-1)+2)
1620:       TEMP_AR(3)= Q(3*(at1-1)+3);TEMP_AR(4)= Q(3*(at2-1)+1)1599:       TEMP_AR(3)= Q(3*(at1-1)+3);TEMP_AR(4)= Q(3*(at2-1)+1)
1621:       TEMP_AR(5)= Q(3*(at2-1)+2);TEMP_AR(6)= Q(3*(at2-1)+3)1600:       TEMP_AR(5)= Q(3*(at2-1)+2);TEMP_AR(6)= Q(3*(at2-1)+3)
1622:       TEMP_AR(7)= Q(3*(at3-1)+1);TEMP_AR(8)= Q(3*(at3-1)+2)1601:       TEMP_AR(7)= Q(3*(at3-1)+1);TEMP_AR(8)= Q(3*(at3-1)+2)
1623:       TEMP_AR(9)= Q(3*(at3-1)+3)1602:       TEMP_AR(9)= Q(3*(at3-1)+3)
1624:       Q = 9999.00001603:       Q = 9999.0000
1625:       DO i= 1,31604:       DO i= 1,3
1626:         Q(3*(at1-1)+i)= TEMP_AR(i)1605:         Q(3*(at1-1)+i)= TEMP_AR(i)
1627:         Q(3*(at2-1)+i)= TEMP_AR(3+i)1606:         Q(3*(at2-1)+i)= TEMP_AR(3+i)
1628:         Q(3*(at3-1)+i)= TEMP_AR(6+i)1607:         Q(3*(at3-1)+i)= TEMP_AR(6+i)
1629:       ENDDO1608:       ENDDO
1630:       CALL CHARMMBILDC(1, lenic2, Q, BOND1, BOND2,THET1,&1609:       CALL CHARMMBILDC(1, lenic2, Q, BOND1, BOND2,THET1,&
1631:      &     THET2, PHI)1610:      &     THET2, PHI)
1632: !      PRINT*, "interpolated"1611:       DO i =1, NATOMS
1633: !      DO i =1, NATOMS1612:          PRINT '(a4,4f11.5)',ih(m04+i-1),Q(3*(i-1)+1),&
1634: !         PRINT '(a4,4f11.5)',ih(m04+i-1),Q(3*(i-1)+1),&1613:      &          Q(3*(i-1)+2),Q(3*(i-1)+3)
1635: !     &          Q(3*(i-1)+2),Q(3*(i-1)+3)1614:       ENDDO
1636: !      ENDDO 
1637: 1615: 
1638:       END SUBROUTINE TAKESTEPAMDIHED1616:       END SUBROUTINE TAKESTEPAMDIHED
1639: 1617: 
1640: ! ************************************************************************ 
1641:       DOUBLE PRECISION FUNCTION PHIINTERP(ST_PHI, FI_PHI, SFRAC) 
1642:       IMPLICIT NONE 
1643:       DOUBLE PRECISION:: ST_PHI,FI_PHI 
1644:       DOUBLE PRECISION:: PHI 
1645:       DOUBLE PRECISION:: SFRAC  
1646:  
1647:       IF (ST_PHI*FI_PHI .GE. 0.D0) THEN  
1648:          PHI = SFRAC*ST_PHI + (1.D0-SFRAC)*FI_PHI 
1649:       ELSEIF (ST_PHI .GE. 0.D0 .AND. FI_PHI .LT. 0.D0) THEN 
1650:          IF (DABS(ST_PHI-FI_PHI).LE.180.D0) THEN 
1651:             PHI = SFRAC*ST_PHI + (1.D0-SFRAC)*FI_PHI 
1652:          ELSE 
1653:             IF (ST_PHI.GT.DABS(FI_PHI)) THEN 
1654:                ST_PHI=ST_PHI-360.D0 
1655:             ELSE 
1656:                FI_PHI=FI_PHI+360.D0 
1657:             ENDIF 
1658:             PHI = SFRAC*ST_PHI + (1.D0-SFRAC)*FI_PHI 
1659:          ENDIF 
1660:       ELSEIF (ST_PHI .LT. 0.D0 .AND. FI_PHI .GE. 0.D0) THEN 
1661:          IF (DABS(ST_PHI-FI_PHI).LE.180.D0) THEN 
1662:             PHI = SFRAC*ST_PHI + (1.D0-SFRAC)*FI_PHI 
1663:          ELSE 
1664:             IF (DABS(ST_PHI).GT.FI_PHI) THEN 
1665:                ST_PHI=ST_PHI+360.D0 
1666:             ELSE 
1667:                FI_PHI=FI_PHI-360.D0 
1668:             ENDIF 
1669:             PHI = SFRAC*ST_PHI + (1.D0-SFRAC)*FI_PHI 
1670:          ENDIF 
1671:       ENDIF !msb50 - take next two lines out of if 
1672:       IF (PHI.GT.180.D0) PHI=PHI-360.D0 
1673:       IF (PHI.LE.-180.D0) PHI=PHI+360.D0 
1674:       PHIINTERP = PHI 
1675:       RETURN 
1676:       END FUNCTION PHIINTERP 
1677: 1618: 
1678: ! *************************************************************************1619: ! *************************************************************************
1679: 1620: 
1680: 1621: 
1681:       SUBROUTINE GETSEEDATM(n1,n2,n3,at1,at2,at3, resnum)1622:       SUBROUTINE GETSEEDATM(n1,n2,n3,at1,at2,at3, resnum)
1682:       use modamber91623:       use modamber9
1683:       use commons, only : natoms1624:       use commons, only : natoms
1684:       INTEGER :: i1625:       INTEGER :: i
1685:       INTEGER :: atomOffset, nat_res,k1626:       INTEGER :: atomOffset, nat_res,k
1686:       INTEGER, INTENT(OUT)::at1, at2,at31627:       INTEGER, INTENT(OUT)::at1, at2,at3
1696:       IF (resnum .GT. nres) THEN1637:       IF (resnum .GT. nres) THEN
1697:          PRINT*, "ERROR IN GETSEEDATM - max no of residues", nres1638:          PRINT*, "ERROR IN GETSEEDATM - max no of residues", nres
1698:          STOP1639:          STOP
1699:       ENDIF1640:       ENDIF
1700: 1641: 
1701:       IF (resnum == 0) THEN1642:       IF (resnum == 0) THEN
1702:           i = 11643:           i = 1
1703:       ELSE1644:       ELSE
1704:           i = resnum1645:           i = resnum
1705:           atomOffset = ix(i02+i-1)-11646:           atomOffset = ix(i02+i-1)-1
 1647:           print*, "atomoffset", atomOffset
1706:       ENDIF 1648:       ENDIF 
1707: 1649: 
1708: 1650: 
1709:       DO WHILE (found .EQV. .FALSE.)1651:       DO WHILE (found .EQV. .FALSE.)
1710:           count = count +11652:           count = count +1
1711:           IF (count .GT. nres) EXIT1653:           IF (count .GT. nres) EXIT
1712: 1654: 
1713:           nat_res =ix(i02+i)- ix(i02+i-1)1655:           nat_res =ix(i02+i)- ix(i02+i-1)
1714:           DO k=atomOffset+1, atomOffset + nat_res1656:           DO k=atomOffset+1, atomOffset + nat_res
1715:               IF (n1 ==ih(m04+k-1)) THEN1657:               IF (n1 ==ih(m04+k-1)) THEN
1726:      &               (found3.EQV..TRUE.)) THEN1668:      &               (found3.EQV..TRUE.)) THEN
1727:                  found = .TRUE.1669:                  found = .TRUE.
1728:                  EXIT1670:                  EXIT
1729:               ENDIF1671:               ENDIF
1730:           ENDDO1672:           ENDDO
1731: 1673: 
1732:           IF (resnum == 0) THEN1674:           IF (resnum == 0) THEN
1733:              i = i+11675:              i = i+1
1734:              atomOffset = atomOffset + nat_res1676:              atomOffset = atomOffset + nat_res
1735:           ENDIF1677:           ENDIF
1736:       1678:       !this is gonna give a segfault if wrong n1, n2,n3 given
1737: !this is gonna give a segfault if wrong n1, n2,n3 given 
1738:       ENDDO1679:       ENDDO
 1680:       PRINT*, "ats", at1, at2, at3
1739:       END SUBROUTINE GETSEEDATM1681:       END SUBROUTINE GETSEEDATM
1740: 1682: 
1741: 1683: 
1742: ! ************************************************************************1684: ! ************************************************************************
1743: 1685: 
1744: 1686: 
1745: 1687: 
1746:       SUBROUTINE CHGETICVAL(y, BOND1,BOND2,ANGLE1, ANGLE2, DIHEDS,DIHED_ONLY) 1688:       SUBROUTINE CHGETICVAL(y, BOND1,BOND2,ANGLE1, ANGLE2, DIHEDS) 
1747:       use modamber91689:       use modamber9
1748:       use commons, only : natoms1690:       use commons, only : natoms
1749: 1691: 
1750: !      INTEGER, INTENT(IN) :: IC_COORDS(nres*95+5*nres,4)1692: !      INTEGER, INTENT(IN) :: IC_COORDS(nres*95+5*nres,4)
1751: !      CHARACTER(LEN=3), INTENT(IN) :: IC_IMPROP(nres*95+5*nres)1693: !      CHARACTER(LEN=3), INTENT(IN) :: IC_IMPROP(nres*95+5*nres)
1752: 1694: 
1753:       double precision :: y(3*natoms)1695:       double precision :: y(3*natoms)
1754:       logical, intent(in) :: DIHED_ONLY 
1755:       DOUBLE PRECISION, INTENT(OUT) :: DIHEDS(nphia+nphih),&1696:       DOUBLE PRECISION, INTENT(OUT) :: DIHEDS(nphia+nphih),&
1756:      &          ANGLE1(ntheth+ntheta+ngper), ANGLE2(ntheth+ntheta*ngper)1697:      &          ANGLE1(ntheth+ntheta+ngper), ANGLE2(ntheth+ntheta*ngper)
1757:       double precision :: bond1(nbonh+nbona+nbper),bond2(nbonh+nbona+nbper)1698:       double precision :: bond1(nbonh+nbona+nbper),bond2(nbonh+nbona+nbper)
1758:       double precision :: dihed1699:       double precision :: dihed
1759:       logical :: dihed_Side(nphia+nphih),angle_Side(ntheth+ntheta+ngper)1700:       logical :: dihed_Side(nphia+nphih),angle_Side(ntheth+ntheta+ngper)
1760:       double precision :: rij, xji, ant , rij2,ant21701:       double precision :: rij, xji, ant , rij2,ant2
1761:       double precision :: vec_ij(3), vec_kj(3)1702:       double precision :: vec_ij(3), vec_kj(3)
1762:       integer :: i,j, i3, j3, k3, l3, count1703:       integer :: i,j, i3, j3, k3, l3, count
1763: !      INTEGER :: lenic1704: !      INTEGER :: lenic
1764:       double precision pt9991705:       double precision pt999
1771: !         PRINT*, IC_COORDS(i,1),IC_COORDS(i,2),IC_COORDS(i,3),IC_COORDS(i,4)1712: !         PRINT*, IC_COORDS(i,1),IC_COORDS(i,2),IC_COORDS(i,3),IC_COORDS(i,4)
1772:          i3 = (IC_COORDS(i,1)-1)*3; j3 = (IC_COORDS(i,2)-1)*31713:          i3 = (IC_COORDS(i,1)-1)*3; j3 = (IC_COORDS(i,2)-1)*3
1773:          k3 = (IC_COORDS(i,3)-1)*3; l3 = (IC_COORDS(i,4)-1)*31714:          k3 = (IC_COORDS(i,3)-1)*3; l3 = (IC_COORDS(i,4)-1)*3
1774:          IF (IC_IMPROP(i) .EQ. "IMP") THEN1715:          IF (IC_IMPROP(i) .EQ. "IMP") THEN
1775:              rij = 0; rij2 = 01716:              rij = 0; rij2 = 0
1776:              DO j = 1,31717:              DO j = 1,3
1777:                  xij = y(i3+j)-y(k3+j)1718:                  xij = y(i3+j)-y(k3+j)
1778:                  rij = rij + xij**21719:                  rij = rij + xij**2
1779:              ENDDO            1720:              ENDDO            
1780:              rij = SQRT(rij)1721:              rij = SQRT(rij)
 1722:              CALL AMBERANGLE(y, i3, k3, j3, ant)
1781:              CALL AMBERDIHEDR(y, NATOMS, i3/3+1, j3/3+1, k3/3+1, l3/3+1,dihed)1723:              CALL AMBERDIHEDR(y, NATOMS, i3/3+1, j3/3+1, k3/3+1, l3/3+1,dihed)
1782:              IF (.NOT. DIHED_ONLY) THEN !otherwise calculate only dihedrals1724:              CALL AMBERANGLE(y, j3,k3,l3,ant2)
1783:                  CALL AMBERANGLE(y, i3, k3, j3, ant)1725:              DO j = 1,3
1784:                  CALL AMBERANGLE(y, j3,k3,l3,ant2)1726:                  xij = y(k3+j)-y(l3+j)
1785:                  DO j = 1,31727:                  rij2 = rij2 + xij**2
1786:                      xij = y(k3+j)-y(l3+j)1728:              ENDDO
1787:                      rij2 = rij2 + xij**21729:              rij2 = SQRT(rij2)
1788:                  ENDDO1730:              bond1(i) = rij; bond2(i) = rij2
1789:                  rij2 = SQRT(rij2)1731:              angle1(i) = ant; angle2(i) = ant2
1790:                  bond1(i) = rij; bond2(i) = rij2 
1791:                  angle1(i) = ant; angle2(i) = ant2 
1792:              ENDIF 
1793:              diheds(i) = dihed1732:              diheds(i) = dihed
1794: !            PRINT '(i3, a4, 4i3,f8.4,3f8.2 f8.4)',i,"IMP",1733: !            PRINT '(i3, a4, 4i3,f8.4,3f8.2 f8.4)',i,"IMP",
1795: !     &        IC_COORDS(i,1),IC_COORDS(i,2),IC_COORDS(i,3),IC_COORDS(i,4),1734: !     &        IC_COORDS(i,1),IC_COORDS(i,2),IC_COORDS(i,3),IC_COORDS(i,4),
1796: !     &        rij, ant, dihed, ant2, rij21735: !     &        rij, ant, dihed, ant2, rij2
1797:          ELSE1736:          ELSE
1798:             rij = 0; rij2 = 0 1737:             rij = 0; rij2 = 0 
1799:             DO j = 1,31738:             DO j = 1,3
1800:                  xij = y(i3+j)-y(j3+j)1739:                  xij = y(i3+j)-y(j3+j)
1801:                  rij = rij + xij**21740:                  rij = rij + xij**2
1802:              ENDDO1741:              ENDDO
1803:              rij = SQRT(rij)1742:              rij = SQRT(rij)
 1743:              CALL AMBERANGLE(y, i3, j3, k3, ant)
1804:              CALL AMBERDIHEDR(y, NATOMS, i3/3+1, j3/3+1, k3/3+1, l3/3+1,dihed)1744:              CALL AMBERDIHEDR(y, NATOMS, i3/3+1, j3/3+1, k3/3+1, l3/3+1,dihed)
1805:              IF (.NOT. DIHED_ONLY) THEN 1745:              CALL AMBERANGLE(y, j3,k3,l3,ant2)
1806:                  CALL AMBERANGLE(y, i3, j3, k3, ant)1746:              DO j = 1,3
1807:                  CALL AMBERANGLE(y, j3,k3,l3,ant2)1747:                  xij = y(k3+j)-y(l3+j)
1808:                  DO j = 1,31748:                  rij2 = rij2 + xij**2
1809:                      xij = y(k3+j)-y(l3+j)1749:              ENDDO
1810:                      rij2 = rij2 + xij**21750:              rij2 = SQRT(rij2)
1811:                  ENDDO1751:              bond1(i) = rij; bond2(i) = rij2
1812:                  rij2 = SQRT(rij2)1752:              angle1(i) = ant; angle2(i) = ant2
1813:                  bond1(i) = rij; bond2(i) = rij2 
1814:                  angle1(i) = ant; angle2(i) = ant2 
1815:              ENDIF 
1816:              diheds(i) = dihed1753:              diheds(i) = dihed
1817: !             PRINT '(i3, a4, 4i3, f8.4, 3f8.2, f8.4)',i, "NOT",1754: !             PRINT '(i3, a4, 4i3, f8.4, 3f8.2, f8.4)',i, "NOT",
1818: !     &        IC_COORDS(i,1),IC_COORDS(i,2),IC_COORDS(i,3),IC_COORDS(i,4),1755: !     &        IC_COORDS(i,1),IC_COORDS(i,2),IC_COORDS(i,3),IC_COORDS(i,4),
1819: !     &        rij, ant, dihed, ant2, rij21756: !     &        rij, ant, dihed, ant2, rij2
1820:          ENDIF1757:          ENDIF
1821: !  ENDDO1758: !  ENDDO
1822:       ENDDO1759:       ENDDO
1823:       END SUBROUTINE CHGETICVAL1760:       END SUBROUTINE CHGETICVAL
1824:       1761:       
1825: 1762: 
1888: 1825: 
1889: !**************1826: !**************
1890:       SUBROUTINE GETICCOORDS()1827:       SUBROUTINE GETICCOORDS()
1891: 1828: 
1892:       use modamber91829:       use modamber9
1893:       use commons, only : NATOMS1830:       use commons, only : NATOMS
1894:       INTEGER :: i, j , a_index1831:       INTEGER :: i, j , a_index
1895:       INTEGER :: atomOffset,  nat_res, var1832:       INTEGER :: atomOffset,  nat_res, var
1896: !      INTEGER :: IC_COORDS(nres*95+5*nres,4)1833: !      INTEGER :: IC_COORDS(nres*95+5*nres,4)
1897:       INTEGER :: IC_START(23) !for each amino acid:where info in ICATOMS starts1834:       INTEGER :: IC_START(23) !for each amino acid:where info in ICATOMS starts
1898:       INTEGER :: IC_STARN(3), IC_STARC(3), IC_STARNPRO(3) !dummies for N/C terminals1835:       INTEGER :: IC_STARN(3), IC_STARC(3) !dummies for N/C terminals
1899:       CHARACTER(LEN=4) :: IC_ATOMS(1250)!(1121) !IMP if improper1836:       CHARACTER(LEN=4) :: IC_ATOMS(1250)!(1121) !IMP if improper
1900:       CHARACTER(LEN=4) :: IC_NTERM(15), IC_CTERM(10), IC_NPRO(10)1837:       CHARACTER(LEN=4) :: IC_NTERM(15), IC_CTERM(10)
1901:       CHARACTER(LEN=3) :: AminAcid(22)1838:       CHARACTER(LEN=3) :: AminAcid(22)
1902: !      CHARACTER(LEN=3) :: IC_IMPROP(nres*95+5*nres)1839: !      CHARACTER(LEN=3) :: IC_IMPROP(nres*95+5*nres)
1903: !      # IC per acid/25,95,45,35, 30, 60, 50,10,60,65,60,70,70,80,60,75,50,1840: !      # IC per acid/25,95,45,35, 30, 60, 50,10,60,65,60,70,70,80,60,75,50,
1904: !                     30, 45,95, 80,551841: !                     30, 45,95, 80,55
1905:       data AminAcid /"ALA","ARG","ASN","ASP","CYS","GLN","GLU","GLY",      &1842:       data AminAcid /"ALA","ARG","ASN","ASP","CYS","GLN","GLU","GLY",      &
1906:      &               "HIE","HIP","HID","ILE","LEU","LYS","MET","PHE",      & 1843:      &               "HIE","HIP","HID","ILE","LEU","LYS","MET","PHE",      & 
1907:      &               "PRO","SER","THR","TRP","TYR","VAL"/                   1844:      &               "PRO","SER","THR","TRP","TYR","VAL"/                   
1908:       data IC_START /0,25,120,165,200,230,290,340,350,410,475, 535, &1845:       data IC_START /0,25,120,165,200,230,290,340,350,410,475, 535, &
1909:      &                605,675,760, 820, 895,945,975,1020,1115,1195,1250/1846:      &                605,675,760, 820, 895,945,975,1020,1115,1195,1250/
1910:       data IC_STARN /1,0,3 /, IC_STARC /1,0,2/, IC_STARNPRO /1,0,2/1847:       data IC_STARN /1,0,3 /, IC_STARC /1,0,2/
1911:       data IC_NPRO  /"H2", "CA", "N", "CD", "IMP",                        & 
1912:      &                "H3", "CA", "N", "H2", "IMP"/ 
1913:       data IC_NTERM /"H1","N","CA","C","NOT",&1848:       data IC_NTERM /"H1","N","CA","C","NOT",&
1914:      &               "H2","CA","N","H1","IMP",&1849:      &               "H2","CA","N","H1","IMP",&
1915:      &               "H3","CA","N","H2","IMP"/1850:      &               "H3","CA","N","H2","IMP"/
1916:       data IC_CTERM /"N", "CA", "C","OXT","NOT",&1851:       data IC_CTERM /"N", "CA", "C","OXT","NOT",&
1917:      &               "OXT", "CA","C","O", "IMP"/1852:      &               "OXT", "CA","C","O", "IMP"/
1918:       data IC_ATOMS /"N", "C", "CA", "CB", "IMP",                         &1853:       data IC_ATOMS /"N", "C", "CA", "CB", "IMP",                         &
1919:      &                 "N", "C", "CA", "HA", "IMP",                       &1854:      &                 "N", "C", "CA", "HA", "IMP",                       &
1920:      &                 "C", "CA", "CB", "HB1", "NOT",                     &1855:      &                 "C", "CA", "CB", "HB1", "NOT",                     &
1921:      &                 "HB1", "CA", "CB", "HB2", "IMP",                   &1856:      &                 "HB1", "CA", "CB", "HB2", "IMP",                   &
1922:      &                 "HB1", "CA", "CB", "HB3", "IMP", &!end ala          1857:      &                 "HB1", "CA", "CB", "HB3", "IMP", &!end ala          
2193: !              PRINT*, ih(m02+i-1), "A_index", a_index2128: !              PRINT*, ih(m02+i-1), "A_index", a_index
2194:               CALL IC_NO2AT(a_index,nat_res,nat_res2,atomOffset,              &2129:               CALL IC_NO2AT(a_index,nat_res,nat_res2,atomOffset,              &
2195:      &                     IC_START,IC_ATOMS)2130:      &                     IC_START,IC_ATOMS)
2196: 2131: 
2197: !              PRINT*, IC_COORDS(lenic-1,1),IC_COORDS(lenic-1,2),2132: !              PRINT*, IC_COORDS(lenic-1,1),IC_COORDS(lenic-1,2),
2198: !     &           IC_COORDS(lenic-1,3),IC_COORDS(lenic-1,4)2133: !     &           IC_COORDS(lenic-1,3),IC_COORDS(lenic-1,4)
2199:               EXIT2134:               EXIT
2200: 2135: 
2201:           ELSEIF (ih(m02+i-1) .EQ. ("N"//AminAcid(j))) THEN !NTERM2136:           ELSEIF (ih(m02+i-1) .EQ. ("N"//AminAcid(j))) THEN !NTERM
2202:               a_index=j2137:               a_index=j
2203:               IF (AminAcid(j) .EQ. "PRO") THEN !proline is different2138:               CALL IC_NO2AT(a_index,nat_res,0 ,atomOffset,                   &
2204:                  CALL IC_NO2AT(a_index,nat_res,0 ,atomOffset,IC_STARNPRO,IC_NPRO) 
2205:               ELSE 
2206:                  CALL IC_NO2AT(a_index,nat_res,0 ,atomOffset,                   & 
2207:      &                     IC_STARN,IC_NTERM)2139:      &                     IC_STARN,IC_NTERM)
2208:               ENDIF 
2209: !              PRINT*, IC_COORDS(1,1),IC_COORDS(1,2),IC_COORDS(1,3),IC_COORDS(1,4)2140: !              PRINT*, IC_COORDS(1,1),IC_COORDS(1,2),IC_COORDS(1,3),IC_COORDS(1,4)
2210: !              PRINT*, ih(m02+i-1), "A_index", a_index2141: !              PRINT*, ih(m02+i-1), "A_index", a_index
2211:               CALL IC_NO2AT(a_index,nat_res, nat_res2,atomOffset,            &2142:               CALL IC_NO2AT(a_index,nat_res, nat_res2,atomOffset,            &
2212:      &                     IC_START,IC_ATOMS)2143:      &                     IC_START,IC_ATOMS)
2213:               EXIT2144:               EXIT
2214:           ELSEIF (ih(m02+i-1) .EQ. ("C"//AminAcid(j))) THEN !CTERM2145:           ELSEIF (ih(m02+i-1) .EQ. ("C"//AminAcid(j))) THEN !CTERM
2215:               a_index =j2146:               a_index =j
2216: !              PRINT*, ih(m02+i-1), "A_index", a_index2147: !              PRINT*, ih(m02+i-1), "A_index", a_index
2217:               CALL IC_NO2AT(a_index,nat_res, 0 ,atomOffset,                  &2148:               CALL IC_NO2AT(a_index,nat_res, 0 ,atomOffset,                  &
2218:      &                     IC_START,IC_ATOMS)2149:      &                     IC_START,IC_ATOMS)


r14089/bhinterp.f90 2017-01-21 10:32:30.831188235 +0000 r14088/bhinterp.f90 2017-01-21 10:32:35.271188235 +0000
 20: !  basin-hopping global optimisation based on the actual PE plus the 20: !  basin-hopping global optimisation based on the actual PE plus the
 21: !  energy of two srings connected to CSTART and CFINISH 21: !  energy of two srings connected to CSTART and CFINISH
 22: !  Assume we only use the LBFGS minimiser. 22: !  Assume we only use the LBFGS minimiser.
 23: !  Allow a different force constant, BHK, from NEB phase. 23: !  Allow a different force constant, BHK, from NEB phase.
 24: ! 24: !
 25: SUBROUTINE BHINTERP(CSTART,CFINISH,NCOORDS,NATOMS,INTERPCOORDS,SUCCESS,DSTART,DFINISH,FINALENERGY,ESTART,EFINISH,DENDPT) 25: SUBROUTINE BHINTERP(CSTART,CFINISH,NCOORDS,NATOMS,INTERPCOORDS,SUCCESS,DSTART,DFINISH,FINALENERGY,ESTART,EFINISH,DENDPT)
 26: USE KEY,ONLY : BHSTEPS, BHCONV, BHTEMP, BHSTEPSIZE, BHACCREJ, DEBUG, MUPDATE, BULKT, TWOD, RIGIDBODY, AMBER, & 26: USE KEY,ONLY : BHSTEPS, BHCONV, BHTEMP, BHSTEPSIZE, BHACCREJ, DEBUG, MUPDATE, BULKT, TWOD, RIGIDBODY, AMBER, &
 27:   &            BHK, BFGSSTEPS, BHDEBUG, BHMAXENERGY, GEOMDIFFTOL, EDIFFTOL, PERMDIST, GMAX, BHSFRAC, AMBERT, NABT, & 27:   &            BHK, BFGSSTEPS, BHDEBUG, BHMAXENERGY, GEOMDIFFTOL, EDIFFTOL, PERMDIST, GMAX, BHSFRAC, AMBERT, NABT, &
 28:   &            BHINTERPUSELOWEST, BHCHECKENERGYT, BHSTEPSMIN 28:   &            BHINTERPUSELOWEST, BHCHECKENERGYT, BHSTEPSMIN
 29: USE MODCHARMM,ONLY : CHRMMT, ICINTERPT, CHECKOMEGAT,CHECKCHIRALT,INTMINT,NOCISTRANS,MINOMEGA 29: USE MODCHARMM,ONLY : CHRMMT, ICINTERPT, CHECKOMEGAT,CHECKCHIRALT,INTMINT,NOCISTRANS,MINOMEGA
 30: USE MODAMBER9, ONLY : NOCISTRANSRNA, GOODSTRUCTURE1,GOODSTRUCTURE2,CISARRAY1,& 30: USE MODAMBER9, ONLY : NOCISTRANSRNA, GOODSTRUCTURE1,GOODSTRUCTURE2,CISARRAY1,CISARRAY2, lenic
 31: CISARRAY2, lenic, NICTOT, AMBERICT, AMBSTEPT, nphia, nphih, AMBPERTT 
 32: USE PORFUNCS 31: USE PORFUNCS
 33: USE COMMONS,ONLY: PARAM1,PARAM2,PARAM3,ZSYM 32: USE COMMONS,ONLY: PARAM1,PARAM2,PARAM3,ZSYM
 34: USE MODEFOL 33: USE MODEFOL
 35: USE SPFUNCTS, ONLY : DUMPCOORDS 34: USE SPFUNCTS, ONLY : DUMPCOORDS
 36: IMPLICIT NONE 35: IMPLICIT NONE
 37: INTEGER, INTENT(IN) :: NCOORDS 36: INTEGER, INTENT(IN) :: NCOORDS
 38: DOUBLE PRECISION CSTART(NCOORDS), CFINISH(NCOORDS) 37: DOUBLE PRECISION CSTART(NCOORDS), CFINISH(NCOORDS)
 39: DOUBLE PRECISION COORDS(NCOORDS), VNEW(NCOORDS), COORDSO(NCOORDS), DSTART, DFINISH, RMAT(3,3), ESTRING, ESPREV, RANDOM, DPRAND, & 38: DOUBLE PRECISION COORDS(NCOORDS), VNEW(NCOORDS), COORDSO(NCOORDS), DSTART, DFINISH, RMAT(3,3), ESTRING, ESPREV, RANDOM, DPRAND, &
 40:   &              BESTCOORDS(NCOORDS), EREAL, RMS2, P0, ETOTAL, EBEST, EPREV, RMS, ENERGY, DELTAX(NCOORDS) 39:   &              BESTCOORDS(NCOORDS), EREAL, RMS2, P0, ETOTAL, EBEST, EPREV, RMS, ENERGY, DELTAX(NCOORDS)
 41: DOUBLE PRECISION INTERPCOORDS(NCOORDS), FINALENERGY, SFRAC, DSTARTSAVE, DFINISHSAVE, ESTART, EFINISH, SLIMIT, FLIMIT 40: DOUBLE PRECISION INTERPCOORDS(NCOORDS), FINALENERGY, SFRAC, DSTARTSAVE, DFINISHSAVE, ESTART, EFINISH, SLIMIT, FLIMIT
 42: DOUBLE PRECISION DIST, DIST2, GMAXSAVE, DENDPT 41: DOUBLE PRECISION DIST, DIST2, GMAXSAVE, DENDPT
 43: DOUBLE PRECISION,DIMENSION(:), ALLOCATABLE:: PHI_START, PHI_FINISH 
 44: INTEGER ITDONE, NSUCCESS, NFAIL, NSUCCESST, NFAILT, NATOMS, J1, J2, ISTAT, K, NSTEPS, NBISECT, MAXBISECT, NTRY, NCOUNT 42: INTEGER ITDONE, NSUCCESS, NFAIL, NSUCCESST, NFAILT, NATOMS, J1, J2, ISTAT, K, NSTEPS, NBISECT, MAXBISECT, NTRY, NCOUNT
 45: LOGICAL MFLAG, ATEST, SUCCESS, PTEST, HITS, HITF, CHIRALFAIL, AMIDEFAIL, EXITBHRUN 43: LOGICAL MFLAG, ATEST, SUCCESS, PTEST, HITS, HITF, CHIRALFAIL, AMIDEFAIL, EXITBHRUN
 46: LOGICAL KNOWE, KNOWG, KNOWH 44: LOGICAL KNOWE, KNOWG, KNOWH
 47: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 45: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 48:  46: 
 49: !  47: ! 
 50: ! Calculate initial energy 48: ! Calculate initial energy
 51: ! 49: !
 52: IF (CHRMMT.AND.INTMINT) THEN 50: IF (CHRMMT.AND.INTMINT) THEN
 53:    PRINT '(A)','bhinterp> CHARMM internal coodinate minimisation not allowed with BHINTERP' 51:    PRINT '(A)','bhinterp> CHARMM internal coodinate minimisation not allowed with BHINTERP'
 75:   &       -PARAM1*NINT((CFINISH(3*(K-1)+1) - CSTART(3*(K-1)+1))/PARAM1) 73:   &       -PARAM1*NINT((CFINISH(3*(K-1)+1) - CSTART(3*(K-1)+1))/PARAM1)
 76:       DELTAX(3*(K-1)+2)=CFINISH(3*(K-1)+2) - CSTART(3*(K-1)+2) & 74:       DELTAX(3*(K-1)+2)=CFINISH(3*(K-1)+2) - CSTART(3*(K-1)+2) &
 77:   &       -PARAM2*NINT((CFINISH(3*(K-1)+2) - CSTART(3*(K-1)+2))/PARAM2) 75:   &       -PARAM2*NINT((CFINISH(3*(K-1)+2) - CSTART(3*(K-1)+2))/PARAM2)
 78:       DELTAX(3*(K-1)+3)=CFINISH(3*(K-1)+3) - CSTART(3*(K-1)+3) & 76:       DELTAX(3*(K-1)+3)=CFINISH(3*(K-1)+3) - CSTART(3*(K-1)+3) &
 79:   &       -PARAM3*NINT((CFINISH(3*(K-1)+3) - CSTART(3*(K-1)+3))/PARAM3) 77:   &       -PARAM3*NINT((CFINISH(3*(K-1)+3) - CSTART(3*(K-1)+3))/PARAM3)
 80:    ENDDO 78:    ENDDO
 81:    COORDS(1:NCOORDS)=SFRAC*CSTART(1:NCOORDS)+(1.0D0-SFRAC)*DELTAX(1:NCOORDS) 79:    COORDS(1:NCOORDS)=SFRAC*CSTART(1:NCOORDS)+(1.0D0-SFRAC)*DELTAX(1:NCOORDS)
 82: ELSE 80: ELSE
 83:    COORDS(1:NCOORDS)=SFRAC*CSTART(1:NCOORDS)+(1.0D0-SFRAC)*CFINISH(1:NCOORDS) 81:    COORDS(1:NCOORDS)=SFRAC*CSTART(1:NCOORDS)+(1.0D0-SFRAC)*CFINISH(1:NCOORDS)
 84:    IF(CHRMMT.AND.ICINTERPT) CALL ICINTERPOL(COORDS,CSTART,CFINISH,SFRAC) 82:    IF(CHRMMT.AND.ICINTERPT) CALL ICINTERPOL(COORDS,CSTART,CFINISH,SFRAC)
 85:    IF (AMBERT.AND.AMBERICT) THEN 83:    IF (AMBERT.AND.ICINTERPT) THEN
 86: !   IF (AMBERT.AND.ICINTERPT) THEN 
 87:       PRINT*, "AMBER internal coordinate interpolation" 84:       PRINT*, "AMBER internal coordinate interpolation"
 88:       IF (.NOT.ALLOCATED(NICTOT)) THEN 
 89:           CALL SETDIHEAM() 
 90:       ENDIF 
 91:       CALL TAKESTEPAMDIHED(COORDS, CSTART, CFINISH,SFRAC) !msb50 85:       CALL TAKESTEPAMDIHED(COORDS, CSTART, CFINISH,SFRAC) !msb50
 92:    ENDIF 86:    ENDIF
 93: ENDIF 87: ENDIF
 94: IF(BHDEBUG) CALL DUMPCOORDS(COORDS, 'midpoint.xyz', .FALSE.) 88: IF(BHDEBUG) CALL DUMPCOORDS(COORDS, 'midpoint.xyz', .FALSE.)
 95: !bs360 89: !bs360
 96: IF (CHRMMT) THEN 90: IF (CHRMMT) THEN
 97:    AMIDEFAIL=.FALSE. 91:    AMIDEFAIL=.FALSE.
 98:    IF (CHECKOMEGAT) THEN 92:    IF (CHECKOMEGAT) THEN
 99:       CALL CHECKOMEGA(COORDS,AMIDEFAIL) 93:       CALL CHECKOMEGA(COORDS,AMIDEFAIL)
100:       IF(BHDEBUG) PRINT '(A,L5)','           interpolated geometry: AMIDEFAIL  = ',AMIDEFAIL 94:       IF(BHDEBUG) PRINT '(A,L5)','           interpolated geometry: AMIDEFAIL  = ',AMIDEFAIL
142:          PRINT *,'bhinterp> could not find acceptable interpolated minimum'136:          PRINT *,'bhinterp> could not find acceptable interpolated minimum'
143:          RETURN137:          RETURN
144:       ENDIF138:       ENDIF
145:       IF(MOD(NTRY,2).EQ.0) NCOUNT=NCOUNT+1139:       IF(MOD(NTRY,2).EQ.0) NCOUNT=NCOUNT+1
146:       PRINT '(A,F8.3)','bhinterp> retry interpolation, SFRAC= ',SFRAC140:       PRINT '(A,F8.3)','bhinterp> retry interpolation, SFRAC= ',SFRAC
147:       GOTO 10141:       GOTO 10
148:    ENDIF142:    ENDIF
149: ENDIF143: ENDIF
150: 144: 
151: IF (PERMDIST) THEN145: IF (PERMDIST) THEN
152:    !msb50 changed order CSTART, COORDS because this is only for distance so shouldn't matter146:    CALL MINPERMDIST(COORDS,CSTART,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
153:    CALL MINPERMDIST(CSTART,COORDS,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT) 
154:    DSTART=SQRT(DIST)147:    DSTART=SQRT(DIST)
155:    CALL MINPERMDIST(CFINISH,COORDS,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)148:    CALL MINPERMDIST(COORDS,CFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
156:    DFINISH=SQRT(DIST)149:    DFINISH=SQRT(DIST)
157: ELSE150: ELSE
158:    CALL NEWMINDIST(COORDS,CSTART, NATOMS,DSTART, BULKT,TWOD,'AX   ',.TRUE.,RIGIDBODY,DEBUG,RMAT)151:    CALL NEWMINDIST(COORDS,CSTART, NATOMS,DSTART, BULKT,TWOD,'AX   ',.TRUE.,RIGIDBODY,DEBUG,RMAT)
159:    CALL NEWMINDIST(COORDS,CFINISH,NATOMS,DFINISH,BULKT,TWOD,'AX   ',.TRUE.,RIGIDBODY,DEBUG,RMAT)152:    CALL NEWMINDIST(COORDS,CFINISH,NATOMS,DFINISH,BULKT,TWOD,'AX   ',.TRUE.,RIGIDBODY,DEBUG,RMAT)
160: ENDIF153: ENDIF
161: 154: 
162: ESPREV=0.5D0*BHK*(DSTART**2+DFINISH**2)155: ESPREV=0.5D0*BHK*(DSTART**2+DFINISH**2)
163: ! PRINT '(A,5G18.8,I8)','EREAL,ESPREV,DS,DF,SFRAC,ITDONE=',EREAL,ESPREV,DSTART,DFINISH,SFRAC,ITDONE156: ! PRINT '(A,5G18.8,I8)','EREAL,ESPREV,DS,DF,SFRAC,ITDONE=',EREAL,ESPREV,DSTART,DFINISH,SFRAC,ITDONE
164: 157: 
165: IF (.NOT.MFLAG) THEN158: IF (.NOT.MFLAG) THEN
248:    EBEST=EPREV+ESPREV241:    EBEST=EPREV+ESPREV
249: ENDIF242: ENDIF
250: BESTCOORDS(1:NCOORDS)=COORDS(1:NCOORDS)243: BESTCOORDS(1:NCOORDS)=COORDS(1:NCOORDS)
251: NSUCCESS=0; NFAIL=0; NSUCCESST=0; NFAILT=0244: NSUCCESS=0; NFAIL=0; NSUCCESST=0; NFAILT=0
252: EXITBHRUN=.FALSE.245: EXITBHRUN=.FALSE.
253: 246: 
254: DO J1=1,BHSTEPS247: DO J1=1,BHSTEPS
255: !248: !
256: ! Perturb current geometry.249: ! Perturb current geometry.
257: !250: !
258: !  IF (AMBERT .AND. ICINTERPT) THEN251:    IF (AMBERT .AND. ICINTERPT) THEN
259:    IF (AMBERT .AND. AMBSTEPT) THEN 
260:       IF (.NOT.ALLOCATED(NICTOT)) THEN 
261:           CALL SETDIHEAM() 
262:       ENDIF 
263:       IF (AMBPERTT .AND. .NOT.AMBERICT) THEN 
264:    !msb50 - note - Hydrogens may have been swapped around so angles don't  
265:    ! conform to the ones from internal interpolation in all cases (swapped) 
266:  
267:           IF (.NOT.ALLOCATED(PHI_START)) ALLOCATE(PHI_START(nphia +nphih)) 
268:           IF (.NOT.ALLOCATED(PHI_FINISH)) ALLOCATE(PHI_FINISH(nphia+nphih)) 
269:           CALL CHGETICVAL(CSTART, PHI_START, PHI_START,PHI_START,PHI_START,PHI_START, .TRUE.) 
270:           CALL CHGETICVAL(CFINISH, PHI_FINISH, PHI_FINISH,PHI_FINISH,PHI_FINISH,PHI_FINISH, .TRUE.) 
271:           !DO K = 1, 104 
272:           !   PRINT*,k,"start", PHI_START(k),"FINISH",PHI_FINISH(k) 
273:           !ENDDO 
274:           CALL GET_TWISTABLE(PHI_START, PHI_FINISH) 
275:       ENDIF 
276:      PRINT*, "AMBER TAKESTEP"252:      PRINT*, "AMBER TAKESTEP"
 253: !      STOP
 254:      IF (J1.EQ.1) THEN
 255:         CALL SETDIHEAM()
 256:      ENDIF
277:      CALL TAKESTEPAMM(COORDS)257:      CALL TAKESTEPAMM(COORDS)
278:    ELSE IF (CHRMMT) THEN258:    ELSE IF (CHRMMT) THEN
279:       CALL TAKESTEPCH(COORDS) ! changed back due to problems with CYPA!259:       CALL TAKESTEPCH(COORDS) ! changed back due to problems with CYPA!
280: !     DO J2=1,NCOORDS260: !     DO J2=1,NCOORDS
281: !        COORDS(J2)=COORDS(J2)+BHSTEPSIZE*(DPRAND()-0.5D0)*2.0D0261: !        COORDS(J2)=COORDS(J2)+BHSTEPSIZE*(DPRAND()-0.5D0)*2.0D0
282: !     ENDDO262: !     ENDDO
283:    ELSE263:    ELSE
284:       DO J2=1,NCOORDS264:       DO J2=1,NCOORDS
285:          COORDS(J2)=COORDS(J2)+BHSTEPSIZE*(DPRAND()-0.5D0)*2.0D0265:          COORDS(J2)=COORDS(J2)+BHSTEPSIZE*(DPRAND()-0.5D0)*2.0D0
286:       ENDDO266:       ENDDO


r14089/charmmBildc.f 2017-01-21 10:32:31.083188235 +0000 r14088/charmmBildc.f 2017-01-21 10:32:35.583188235 +0000
266:       WC=R*SNT*SNP266:       WC=R*SNT*SNP
267: C267: C
268:       Q(3*(L-1)+1)=Q(3*(K-1)+1)+WA*XA+WB*XB+WC*XC268:       Q(3*(L-1)+1)=Q(3*(K-1)+1)+WA*XA+WB*XB+WC*XC
269:       Q(3*(L-1)+2)=Q(3*(K-1)+2)+WA*YA+WB*YB+WC*YC269:       Q(3*(L-1)+2)=Q(3*(K-1)+2)+WA*YA+WB*YB+WC*YC
270:       Q(3*(L-1)+3)=Q(3*(K-1)+3)+WA*ZA+WB*ZB+WC*ZC270:       Q(3*(L-1)+3)=Q(3*(K-1)+3)+WA*ZA+WB*ZB+WC*ZC
271: C271: C
272:       AT_OK=.TRUE.272:       AT_OK=.TRUE.
273:       RETURN273:       RETURN
274:       END SUBROUTINE CH_CARTCV274:       END SUBROUTINE CH_CARTCV
275: 275: 
 276: C      LOGICAL FUNCTION CHINITIA(ATOM,X,Y,Z)
 277: C      SUBROUTINE CHARMMGETICV(I,J,K,L,T,RIJ,TIJK,PIJKL,TJKL,RKL,
 278: C     1    X,Y,Z)
276: 279: 
277: 280: 
278: C**************************************************************281: C**************************************************************
279: 282: 
280: 283: 
281:       SUBROUTINE CH_SEED(I,J,K,Q,BOND1,BOND2, ANGLE1,284:       SUBROUTINE CH_SEED(I,J,K,Q,BOND1,BOND2, ANGLE1,
282:      &   ANGLE2, PHI)285:      &   ANGLE2, PHI)
283:       use modamber9286:       use modamber9
284:       use commons, only : NATOMS287:       use commons, only : NATOMS
285:       IMPLICIT NONE288:       IMPLICIT NONE
388: C      USE KEY, ONLY: CHRIGIDT391: C      USE KEY, ONLY: CHRIGIDT
389:       IMPLICIT NONE392:       IMPLICIT NONE
390:       DOUBLE PRECISION    Q(3*NATOM)393:       DOUBLE PRECISION    Q(3*NATOM)
391:       DOUBLE PRECISION    CHPMIN,CHPMAX394:       DOUBLE PRECISION    CHPMIN,CHPMAX
392:       INTEGER             ISEG,CHNMIN,CHNMAX,ISEED395:       INTEGER             ISEG,CHNMIN,CHNMAX,ISEED
393: 396: 
394:       CHPMIN=0.2D0397:       CHPMIN=0.2D0
395:       CHPMAX=0.4D0398:       CHPMAX=0.4D0
396:       CHNMIN=1399:       CHNMIN=1
397:       CHNMAX=0400:       CHNMAX=0
398:       !msb50 correction for only twisting angles with difference > perthresh401:       DO ISEG=1,NSEG
399:       IF (AMBPERTT) THEN402:          CHNMAX=CHNMAX+NPHIPSI(ISEG)+NSIDECHAIN(ISEG)
400:             CHNMAX = NTW_ANGLES403:       ENDDO
401:             !PRINT*, "CHNMAX", CHNMAX 
402:       ELSE 
403:          DO ISEG=1,NSEG 
404:             CHNMAX=CHNMAX+NPHIPSI(ISEG)+NSIDECHAIN(ISEG) 
405:          ENDDO 
406:       ENDIF 
407:       CHNMAX=INT(CHNMAX/2.5D0) 404:       CHNMAX=INT(CHNMAX/2.5D0) 
408:       ISEED=-1405:       ISEED=-1
409: 406: 
410:       CALL PERTDIHAM(Q,CHPMIN,CHPMAX,CHNMIN,CHNMAX,ISEED)407:       CALL PERTDIHAM(Q,CHPMIN,CHPMAX,CHNMIN,CHNMAX,ISEED)
411: C      IF (CHRIGIDT .AND. NSEG.GT.1) THEN408: C      IF (CHRIGIDT .AND. NSEG.GT.1) THEN
412: C         CALL MKRIGIDTRANS(COORDS)409: C         CALL MKRIGIDTRANS(COORDS)
413: C         CALL MKRIGIDROT(COORDS)410: C         CALL MKRIGIDROT(COORDS)
414: C      ENDIF411: C      ENDIF
415: 412: 
416:       RETURN413:       RETURN
433:       INTEGER::           ATOT,A,B,C,I1,J1,IICD, III430:       INTEGER::           ATOT,A,B,C,I1,J1,IICD, III
434:       INTEGER::           RESNUMseg,J3,ISEG,IRES,TOTPHIPSI,TOTSIDECHAIN431:       INTEGER::           RESNUMseg,J3,ISEG,IRES,TOTPHIPSI,TOTSIDECHAIN
435:       INTEGER::           CHNMIN,CHNMAX,ISEED432:       INTEGER::           CHNMIN,CHNMAX,ISEED
436: !     DOUBLE PRECISION::    X(NATOMS),Y(NATOMS),Z(NATMS)433: !     DOUBLE PRECISION::    X(NATOMS),Y(NATOMS),Z(NATMS)
437:       LOGICAL::   TPP(NATOMS),TS(NATOMS),TO(NATOMS),TT434:       LOGICAL::   TPP(NATOMS),TS(NATOMS),TO(NATOMS),TT
438:       INTEGER::   II,JJ,KK,LL, lenic2435:       INTEGER::   II,JJ,KK,LL, lenic2
439:       DOUBLE PRECISION :: BOND1(nbonh+nbona+nbper), BOND2(nbonh+nbona+nbper),436:       DOUBLE PRECISION :: BOND1(nbonh+nbona+nbper), BOND2(nbonh+nbona+nbper),
440:      1       THET1(ntheth+ntheta+ngper), THET2(ntheth+ntheta+ngper),437:      1       THET1(ntheth+ntheta+ngper), THET2(ntheth+ntheta+ngper),
441:      2       PHI(nphia+nphih)438:      2       PHI(nphia+nphih)
442:       DOUBLE PRECISION BB1,BB2,TT1,TT2,PP439:       DOUBLE PRECISION BB1,BB2,TT1,TT2,PP
443:       DOUBLE PRECISION :: ANGLE_SCALE(lenic)440:       INTEGER:: count
444:       INTEGER:: count, P_RESPAR, P_RESMAX 
445:       DOUBLE PRECISION:: ANGMAX,ANGMIN 
446:       DOUBLE PRECISION :: PROBABIL 
447:       DOUBLE PRECISION :: slope 
448: 441: 
449:       count=0 
450: C      write(*,*)'CHPMIN,CHPMAX,CHNMIN,CHNMAX= ',CHPMIN,CHPMAX,CHNMIN,CHNMAX442: C      write(*,*)'CHPMIN,CHPMAX,CHNMIN,CHNMAX= ',CHPMIN,CHPMAX,CHNMIN,CHNMAX
451:  
452:       !msb50 - preliminary -  how many residues have a lin de/increasing probability 
453:        P_RESPAR = 6 !number of residues for which twisting prob decreases 
454:        P_RESMAX = 15 
455:       !ANGMAX, ANGMIN for angles rescaling - twist more when at chain end 
456:        ANGMAX= 1.2 
457:        ANGMIN= 0.6 
458:  
459: C  fill IC table with actual Cartesians 
460:       CALL CHGETICVAL(Q,BOND1, BOND2, THET1, THET2, PHI, .FALSE.) 
461:       443:       
 444: C  fill IC table with actual Cartesians
 445:       CALL CHGETICVAL(Q,BOND1, BOND2, THET1, THET2, PHI)
462: C initialise random nmber generator with input seed446: C initialise random nmber generator with input seed
463:       IF(ISEED.GT.-1) CALL SDPRND(ISEED)447:       IF(ISEED.GT.-1) CALL SDPRND(ISEED)
464: C448: C
465: C will be sent back to 192 if too many or too few dihedrals are altered449: C will be sent back to 192 if too many or too few dihedrals are altered
466: C as determined by CHNMIN and CHNMAX450: C as determined by CHNMIN and CHNMAX
467: 192   CONTINUE451: 192   CONTINUE
468: 452: 
469:       count = count +1 !too avoid endless loop453:       count = count +1 !too avoid endless loop
470: C454: C
471: C phi/psi angles, omega and sidechain dihedrals are stored in separate lists455: C phi/psi angles, omega and sidechain dihedrals are stored in separate lists
472: C456: C
473: C phi/psi angles457: C phi/psi angles
474:       B=0458:       B=0
475:       ATOT=0459:       ATOT=0
476:       DO ISEG=1,NSEG !not at the moment - check twist.src for amendments460:       DO ISEG=1,NSEG !not at the moment - check twist.src for amendments
477:          DO A=1,NPHIPSI(ISEG)461:          DO A=1,NPHIPSI(ISEG)
478:             TPP(ATOT+A)=.FALSE.462:             TPP(ATOT+A)=.FALSE.
479: C463: C
480: C  Calculate P, the probability of twisting464: C  Calculate P, the probability of twisting
481: C465: C
482:             IF (AMBOLDPERTT) THEN466:             IF (REAL(A).LE.(0.5*NPHIPSI(ISEG))) THEN
483:                 IF (REAL(A).LE.(0.5*(NPHIPSI(ISEG)+1))) THEN467:                P=CHPMAX-A*((CHPMAX-CHPMIN)/(NPHIPSI(ISEG)*0.5))
484:                    !P=CHPMAX-A*((CHPMAX-CHPMIN)/(NPHIPSI(ISEG)*0.5)) 
485:                    slope = (CHPMAX-CHPMIN)/(0.5*(NPHIPSI(ISEG)-1)-1) 
486:                    P=-slope*A+(CHPMAX+slope) 
487:                 ELSE 
488:                    !P=CHPMIN+(A-0.5*NPHIPSI(ISEG))*((CHPMAX-CHPMIN)/(NPHIPSI(ISEG)*0.5)) 
489:                    slope = (CHPMAX-CHPMIN)/(0.5*(NPHIPSI(ISEG)+1)) 
490:                    P = slope*A +(CHPMAX-slope*NPHIPSI(ISEG)) 
491:                 END IF 
492:             ELSE468:             ELSE
493:             !msb50469:                P=CHPMIN+(A-0.5*NPHIPSI(ISEG))*((CHPMAX-CHPMIN)/(NPHIPSI(ISEG)*0.5))
494:             !take 2* P_RESMAX as NPHIPSI(ISEG) = 2*RESNUMseg(ISEG) 470:             END IF
495:             !  one phi, one psi angle per chain 
496:                 P = PROBABIL(A,NPHIPSI(ISEG), CHPMAX,CHPMIN,2*P_RESMAX,2*P_RESPAR) 
497:                !msb50 - calculation of P depends on position of angle in chain only 
498:                ! can't rely on TW_ANGLES being space out equally along the chain --> 
499:                ! would not have P=CHNMIN in the middle of the chain if P= P(NTW_ANGLES) 
500:                 IF (AMBPERTT) THEN 
501:                    !PRINT '(a15, 2f7.3)', "AMBERPERT probs",P,  TW_DIFFP(PHIPSI(A)) 
502:                    P = (P + TW_DIFFP(PHIPSI(A)))/2!rescaled prob+difference prob-then rescale 
503:                 ENDIF 
504:  
505:                 ! Calculate scaling factor for angle 
506:                 ANGLE_SCALE(PHIPSI(ATOT+A)) = PROBABIL(A,NPHIPSI(ISEG), 
507:      &                        ANGMAX,ANGMIN,2*P_RESMAX,2*P_RESPAR) 
508:             ENDIF 
509: C            IF(BHDEBUG) WRITE(*,*)'P phipsi = ',P,ATOT+A471: C            IF(BHDEBUG) WRITE(*,*)'P phipsi = ',P,ATOT+A
510:             MYRANDOM=DPRAND()472:             MYRANDOM=DPRAND()
511:             !PRINT '(i4,a24,2f7.3)',PHIPSI(A), "nphipsi probabilities", MYRANDOM,P 
512:             MYRANDOM=DPRAND() 
513:             IF (MYRANDOM.LT.P) THEN473:             IF (MYRANDOM.LT.P) THEN
514:                TPP(ATOT+A)=.TRUE.474:                TPP(ATOT+A)=.TRUE.
515:                B=B+1475:                B=B+1
516:             ENDIF476:             ENDIF
517:          ENDDO477:          ENDDO
518:          ATOT=ATOT+NPHIPSI(ISEG)478:          ATOT=ATOT+NPHIPSI(ISEG)
519:       ENDDO479:       ENDDO
520:       TOTPHIPSI=ATOT480:       TOTPHIPSI=ATOT
521: C      IF(BHDEBUG) WRITE(*,*)'PHIPSI: TOT=',ATOT481: C      IF(BHDEBUG) WRITE(*,*)'PHIPSI: TOT=',ATOT
522: 482: 
558: C  this means that all dihedrals in the same sidechain have the same P518: C  this means that all dihedrals in the same sidechain have the same P
559: C519: C
560:             IICD=TW_SIDECHAIN(ATOT+A)520:             IICD=TW_SIDECHAIN(ATOT+A)
561:             j3 = IC_COORDS(IICD,2)521:             j3 = IC_COORDS(IICD,2)
562:             DO III = 1,NRES522:             DO III = 1,NRES
563:                IF (ix(i02+III-1).LE.j3 .AND.(ix(i02+III).GT.j3)) THEN523:                IF (ix(i02+III-1).LE.j3 .AND.(ix(i02+III).GT.j3)) THEN
564:                ires = III524:                ires = III
565:                EXIT525:                EXIT
566:             ENDIF526:             ENDIF
567:             ENDDO527:             ENDDO
568:             IF (ISEG .GT. 1) THEN !NICTOT different from CHARMM!!528:             IRES=IRES-NICTOT(ISEG) !no of residue in segment
569:                RESNUMseg=NICTOT(ISEG)-NICTOT(ISEG-1)     ! number of residues in this segment529:             RESNUMseg=NICTOT(ISEG+1)-NICTOT(ISEG)     ! number of residues in this segment
570:                IRES=IRES-NICTOT(ISEG-1) !no of residue in segment 
571:             ELSE  
572:                RESNUMseg=NICTOT(ISEG) 
573:             ENDIF 
574: C            IF(BHDEBUG) WRITE(*,*)'IICD, JAR, IRES, RESNUM : ',IICD,JAR,IRES,RESNUM530: C            IF(BHDEBUG) WRITE(*,*)'IICD, JAR, IRES, RESNUM : ',IICD,JAR,IRES,RESNUM
575:             IF (AMBOLDPERTT) THEN531:             IF (REAL(IRES).LE.(0.5*RESNUMseg)) THEN
576:                 IF (REAL(IRES).LE.(0.5*(RESNUMseg+1))) THEN532:                P=CHPMAX-IRES*(CHPMAX-CHPMIN)/(RESNUMseg*0.5)
577:                    P=CHPMAX-IRES*(CHPMAX-CHPMIN)/(RESNUMseg*0.5) 
578:                    slope = (CHPMAX-CHPMIN)/(0.5*(RESNUMseg-1)-1) !msb50 
579:                    P=-slope*IRES+(CHPMAX+slope) 
580:                 ELSE 
581:                    P=CHPMIN+(IRES-0.5*RESNUMseg)*(CHPMAX-CHPMIN)/(RESNUMseg*0.5) 
582:                    slope = (CHPMAX-CHPMIN)/(0.5*(RESNUMseg+1)) 
583:                    P = slope*IRES +(CHPMAX-slope*RESNUMseg) 
584:                 END IF 
585:             ELSE533:             ELSE
586:                 P = PROBABIL(IRES,RESNUMseg,CHPMAX,CHPMIN,P_RESMAX, P_RESPAR)534:                P=CHPMIN+(IRES-0.5*RESNUMseg)*(CHPMAX-CHPMIN)/(RESNUMseg*0.5)
587:                 IF (AMBPERTT) THEN535:             END IF
588:                    !PRINT '(a20, 2f7.3)', "AMBERPERT probs", P, TW_DIFFP(IICD) 
589:                    P = (P + TW_DIFFP((IICD)))/2 !rescaled prob + difference prob - then rescale 
590:                 ENDIF 
591:                 !  Calculate the angle scaling factor 
592:                 ANGLE_SCALE(IICD) = PROBABIL(IRES,RESNUMseg, 
593:      &                        ANGMAX,ANGMIN,P_RESMAX,P_RESPAR) 
594:             ENDIF 
595:             !PRINT '(i5,a20,i3,a,2f10.7)',IICD, "side probabs for", ires,":", MYRANDOM,P 
596:  
597:             MYRANDOM=DPRAND() 
598:  
599: C            IF(BHDEBUG) WRITE(*,*)'P sidechain =',P,ATOT+A536: C            IF(BHDEBUG) WRITE(*,*)'P sidechain =',P,ATOT+A
600:             MYRANDOM=DPRAND()537:             MYRANDOM=DPRAND()
601:             IF (MYRANDOM.LT.P) THEN538:             IF (MYRANDOM.LT.P) THEN
602:                TS(ATOT+A)=.TRUE.539:                TS(ATOT+A)=.TRUE.
603:                B=B+1540:                B=B+1
604:             ENDIF541:             ENDIF
605:          ENDDO542:          ENDDO
606:          ATOT=ATOT+NSIDECHAIN(ISEG)543:          ATOT=ATOT+NSIDECHAIN(ISEG)
607:       ENDDO544:       ENDDO
608:       TOTSIDECHAIN=ATOT545:       TOTSIDECHAIN=ATOT
609: C      IF(BHDEBUG) WRITE(*,*)'SIDECHAIN: TOT=',ATOT546: C      IF(BHDEBUG) WRITE(*,*)'SIDECHAIN: TOT=',ATOT
610: C547: C
611: C      shifting b dihedrals, should be NTEST1 < B < NTEST2548: C      shifting b dihedrals, should be NTEST1 < B < NTEST2
612:        IF (B.LT.CHNMIN .OR. B.GT.CHNMAX) THEN549:        IF (B.LT.CHNMIN .OR. B.GT.CHNMAX) THEN
613:          WRITE (*,'(A)') 'Too many dihedrals shifted - retrying'550:          WRITE (*,'(A)') 'Too many dihedrals shifted - retrying'
614:          IF (count<11) THEN551:          IF (count<100000) THEN
615:             GOTO 192552:             GOTO 192
616:          ELSEIF (11.LE.count .AND. count.LE.10000) THEN553:          ELSE 
617:             CHNMAX = CHNMAX +1 !msb50 - quite often tries to twist too many      
618:             !if AMBPERTONLY as CHNMAX = NTW_ANGLES << NPHIPSI+NSIDE 
619:             GOTO 192 
620:          ELSE 
621:              PRINT*, "loop in perdiham failed"554:              PRINT*, "loop in perdiham failed"
622:              STOP555:              STOP
623:          ENDIF556:          ENDIF
624:        ENDIF557:        ENDIF
625: 558: 
626: C  twisting phi/psi angles559: C  twisting phi/psi angles
627:        ATOT=0560:        ATOT=0
628:        DO ISEG=1,NSEG561:        DO ISEG=1,NSEG
629:           DO A=1,NPHIPSI(ISEG)562:           DO A=1,NPHIPSI(ISEG)
630:              IF (TPP(ATOT+A)) THEN563:              IF (TPP(ATOT+A)) THEN
631:                 IICD=PHIPSI(ATOT+A)564:                 IICD=PHIPSI(ATOT+A)
632: C                IF(BHDEBUG) PRINT *,'PERTDIHE> changing phipsi ',IICD565: C                IF(BHDEBUG) PRINT *,'PERTDIHE> changing phipsi ',IICD
633:                 IF (AMBOLDPERTT) THEN566:                 ANGLE=(DPRAND()-0.5)*2.0*BHSTEPSIZE
634:                     ANGLE=(DPRAND()-0.5)*2.0*BHSTEPSIZE 
635:                     !PRINT '(a20,i4,2f10.3)', "phipsi change", IICD,PHI(IICD),PHI(IICD)+ANGLE 
636:                 ELSE 
637:                     ANGLE=((DPRAND()-0.5)*2.0*BHSTEPSIZE)*ANGLE_SCALE(IICD) 
638:                     !PRINT '(i4,a6,f7.3,a6,f10.3)',IICD,"SCALE",ANGLE_SCALE(IICD),"ANGLE", ANGLE 
639:                 ENDIF 
640:                 PHI(IICD) = PHI(IICD) + ANGLE567:                 PHI(IICD) = PHI(IICD) + ANGLE
641:              ENDIF568:              ENDIF
642:           ENDDO569:           ENDDO
643:           ATOT=ATOT+A570:           ATOT=ATOT+A
644:        ENDDO571:        ENDDO
645: C572: C
646: C  twisting amide bond`573: C  twisting amide bond
647:        IF (TOMEGAC) THEN574:        IF (TOMEGAC) THEN
648:           ATOT=0575:           ATOT=0
649:           DO ISEG=1,NSEG576:           DO ISEG=1,NSEG
650:              DO A=1,NOMEGAC(ISEG)577:              DO A=1,NOMEGAC(ISEG)
651:                IF (TO(ATOT+A)) THEN578:                IF (TO(ATOT+A)) THEN
652:                   IICD=OMEGAC(ATOT+A)579:                   IICD=OMEGAC(ATOT+A)
653: C                  IF(BHDEBUG) PRINT *,'PERTDIHE> changing omega ',IICD580: C                  IF(BHDEBUG) PRINT *,'PERTDIHE> changing omega ',IICD
654:                   ANGLE=(DPRAND()-0.5)*2.0*BHSTEPSIZE581:                   ANGLE=(DPRAND()-0.5)*2.0*BHSTEPSIZE
655:                   PHI(IICD) = PHI(IICD) +ANGLE582:                   PHI(IICD) = PHI(IICD) +ANGLE
656:                   !PRINT '(a15,i4,2f10.3)', "omega change",IICD, PHI(IICD)-ANGLE,PHI(IICD) 
657:                ENDIF583:                ENDIF
658:              ENDDO584:              ENDDO
659:              ATOT=ATOT+A585:              ATOT=ATOT+A
660:           ENDDO586:           ENDDO
661:        ENDIF587:        ENDIF
662: C588: C
663: C  twisting sidechains589: C  twisting sidechains
664:        ATOT=0590:        ATOT=0
665:        DO ISEG=1,NSEG591:        DO ISEG=1,NSEG
666:           DO A=1,NSIDECHAIN(ISEG)592:           DO A=1,NSIDECHAIN(ISEG)
667:              IF (TS(ATOT+A)) THEN593:              IF (TS(ATOT+A)) THEN
668:                 IICD=TW_SIDECHAIN(ATOT+A)594:                 IICD=TW_SIDECHAIN(ATOT+A)
669: C                IF(BHDEBUG) PRINT *,'PERTDIHE> changing sidechain ',IICD595: C                IF(BHDEBUG) PRINT *,'PERTDIHE> changing sidechain ',IICD
670:                 IF (AMBOLDPERTT) THEN596:                 ANGLE=(DPRAND()-0.5)*2.0*BHSTEPSIZE
671:                    ANGLE=(DPRAND()-0.5)*2.0*BHSTEPSIZE 
672:                    !PRINT '(a20,i4,2f10.3)', "side change", IICD,PHI(IICD),PHI(IICD)+ANGLE 
673:                 ELSE 
674:                    ANGLE=((DPRAND()-0.5)*2.0*BHSTEPSIZE)*ANGLE_SCALE(IICD) 
675:                    !PRINT '(i4,a6,f7.3,a6,f10.3)',IICD,"SCALE",ANGLE_SCALE(IICD),"ANGLE", ANGLE 
676:                 ENDIF 
677:                 PHI(IICD) = PHI(IICD) +ANGLE597:                 PHI(IICD) = PHI(IICD) +ANGLE
678:              ENDIF598:              ENDIF
679:           ENDDO599:           ENDDO
680:           ATOT=ATOT+A600:           ATOT=ATOT+A
681:        ENDDO601:        ENDDO
682: 602: 
683:        IF(PHI(IICD).LT.-180.0) PHI(IICD)=PHI(IICD)+360.0603:        IF(PHI(IICD).LT.-180.0) PHI(IICD)=PHI(IICD)+360.0
684:        IF(PHI(IICD).GT.180.0) PHI(IICD)=PHI(IICD)-360.0604:        IF(PHI(IICD).GT.180.0) PHI(IICD)=PHI(IICD)-360.0
685:  605:  
686:        CALL CHREBUILD(Q, BOND1, BOND2,THET1,606:        CALL CHREBUILD(Q, BOND1, BOND2,THET1,
687:      &        THET2, PHI)607:      &        THET2, PHI)
688: C      DO II =1, NATOMS608: C      DO II =1, NATOMS
689: C         PRINT '(a4,4f11.5)',ih(m04+II-1),Q(3*(II-1)+1),609: C         PRINT '(a4,4f11.5)',ih(m04+II-1),COORDS(3*(II-1)+1),
690: C     &          Q(3*(II-1)+2),Q(3*(II-1)+3)610: C     &          COORDS(3*(II-1)+2),COORDS(3*(II-1)+3)
691: C      ENDDO611: C      ENDDO
692: 612: 
693: 613: 
694:       END SUBROUTINE PERTDIHAM614:       END SUBROUTINE PERTDIHAM
695: 615: 
696: ! ********************************************************************** 
697:  
698:       DOUBLE PRECISION FUNCTION PROBABIL(x, xmax, CHPMAX,CHPMIN,  
699:      & max_x_func, slope_par) 
700:       IMPLICIT NONE 
701:       DOUBLE PRECISION, INTENT(IN):: CHPMAX, CHPMIN 
702:       INTEGER, INTENT(IN) :: x, xmax, max_x_func, slope_par 
703:       DOUBLE PRECISION :: P, slope 
704:       DOUBLE PRECISION :: P_RESPAR 
705:  
706: ! msb50  - this calculates probability for twisting/scaling of how much angles are twisted 
707: !    if xmax < max_x_func                        else: 
708: !             (a certain no of residues) 
709: !      P |                                         P | 
710: !        |  x           x                            |  x                         x 
711: !        |    x       x                              |    x                     x 
712: !        |      x   x                                |      x                 x        
713: !        |        x                                  |        x x x x x x x x 
714: !        |__________________                         |_____________________________ 
715: !           |     |     |                               |     |             |     | 
716: !           1          xmax                             1                        xmax 
717:  
718:       IF (xmax .GT. max_x_func) THEN 
719:          slope = (CHPMAX-CHPMIN)/(slope_par) 
720:          IF (x.LE.slope_par) THEN 
721:             P = -slope*x+(CHPMAX+slope) 
722:          ELSEIF ((slope_par.LT.x).AND.(x.LE.xmax-slope_par)) THEN 
723:             P = CHPMIN 
724:          ELSE  
725:             P = slope*x +(CHPMAX-slope*xmax) 
726:          ENDIF 
727:       ELSE 
728:             slope = (CHPMAX-CHPMIN)/(0.5*(xmax-1)) 
729:          IF (REAL(x).LE.(0.5*(xmax+1))) THEN 
730:             !P=CHPMAX-x*(CHPMAX-CHPMIN)/(xmax*0.5) 
731:             P=-slope*x+(CHPMAX+slope) 
732:          ELSE 
733:             !P=CHPMIN+(x-0.5*xmax)*(CHPMAX-CHPMIN)/(xmax*0.5) 
734:             P = slope*x +(CHPMAX-slope*xmax) 
735:          ENDIF 
736:       ENDIF 
737:  
738:       PROBABIL = P 
739:       RETURN 
740:       END FUNCTION PROBABIL 
741:  
742: ! ********************************************************************** 
743:  
744: 616: 
745:       SUBROUTINE SETDIHEAM()617:       SUBROUTINE SETDIHEAM()
746:       use modamber9618:       use modamber9
747:       use commons619:       use commons
748:       IMPLICIT NONE620:       IMPLICIT NONE
749:       INTEGER IICD,ISLCT,OLDUSD,I3, J3,K3, L3,IRES,JRES,ISEG,A, NSEGLAST621:       INTEGER IICD,ISLCT,OLDUSD,I3, J3,K3, L3,IRES,JRES,ISEG,A, NSEGLAST
750:       INTEGER AP,AO,AS,AC622:       INTEGER AP,AO,AS,AC
751:       INTEGER i,j623:       INTEGER i,j
752:       INTEGER SEG_START(NATOMS)624:       INTEGER SEG_START(NATOMS)
753:       LOGICAL LPHIPSI,LOMEGAC,LSIDECHAIN,LCHIRAL,LASTSIDECHAIN625:       LOGICAL LPHIPSI,LOMEGAC,LSIDECHAIN,LCHIRAL,LASTSIDECHAIN
754:       INTEGER NPHIPSITOT, NOMEGACTOT,NSIDECHAINTOT, NCHIRALTOT626:       INTEGER NPHIPSITOT, NOMEGACTOT,NSIDECHAINTOT, NCHIRALTOT
755:       CHARACTER*4 TYPEI,TYPEJ,TYPEK627:       CHARACTER*4 TYPEI,TYPEJ,TYPEK
756:       CHARACTER*8 RESLAB628:       CHARACTER*8 RESLAB
757: 629: 
758: ! msb50 - Warning: NICTOT defined differently from CHARMM! 
759: !         CHARMM - residues per segment: NICTOT(ISEG+1) - NICTOT(ISEG)  
760: !              i.e. NICTOT(ISEG) is at which residue new segment starts 
761: !         AMBER  - residues per segment: NICTOT(ISEG) - NICTOT(ISEG-1)  
762: !              i.e. NICTOT(ISEG) is at which residue segment finishes 
763:  
764:       IF (.NOT. ALLOCATED(NICTOT)) ALLOCATE(NICTOT(NATOMS))630:       IF (.NOT. ALLOCATED(NICTOT)) ALLOCATE(NICTOT(NATOMS))
765: 631: 
766:       IF (.NOT. ALLOCATED(IC_COORDS)) THEN632:       IF (.NOT. ALLOCATED(IC_COORDS)) THEN
767:          CALL GETICCOORDS()633:          CALL GETICCOORDS()
768:          PRINT*, "lenic", lenic634:          PRINT*, "lenic", lenic
769:       ENDIF635:       ENDIF
770: 636: 
771:       IF (.NOT. ALLOCATED(IS_SIDECHAIN)) THEN 
772:           ALLOCATE(IS_SIDECHAIN(lenic)) 
773:           DO IICD=1, lenic 
774:              i3 = IC_COORDS(IICD,1); j3 = IC_COORDS(IICD,2) 
775:              k3 = IC_COORDS(IICD,3); l3 = IC_COORDS(IICD,4) 
776:              CALL CHECK_SIDECHAIN(i3,j3,k3,l3,IICD,IS_SIDECHAIN) 
777:           ENDDO 
778:       ENDIF 
779:  
780:       nseg = 0637:       nseg = 0
781:       i = 1638:       i = 1
782:       !SEG_START _ first atom of new segment639:       !SEG_START _ first atom of new segment
783:       SEG_START(1) =1640:       SEG_START(1) =1
784:       DO WHILE (i .LE. NATOMS)641:       DO WHILE (i .LE. NATOMS)
785:       IF (ih(m04+i-1).EQ.'OXT ') THEN !CTerm642:       IF (ih(m04+i-1).EQ.'OXT ') THEN !CTerm
786:          nseg = nseg +1643:          nseg = nseg +1
787:          DO j = 1, nres  !ix(i02+1) = no of atoms in 1st residue +1644:          DO j = 1, nres  !ix(i02+1) 0 no of atoms in 1st residue +1
788:             IF (ix(i02+j-1).LT.i .AND. ix(i02+j).GT.i) THEN645:             IF (ix(i02+j-1).LT.i .AND. ix(i02+j).GT.i) THEN
789:             NICTOT(nseg) = j !number of residues prior to segment + in segment646:             NICTOT(nseg) = j !number of residues per segment
790:             EXIT647:             EXIT
791:             ENDIF648:             ENDIF
792:          ENDDO649:          ENDDO
793: 650: 
794:          i=i+1651:          i=i+1
795:          IF (i .LT. NATOMS) THEN652:          IF (i .LT. NATOMS) THEN
796:             IF (ih(m04+i-1).EQ.'N') THEN653:             IF (ih(m04+i-1).EQ.'N') THEN
797:                SEG_START(nseg+1) = i654:                SEG_START(nseg+1) = i
798:             ELSEIF (ih(m04+i).EQ.'CH3 ') THEN !starts at i+1+1655:             ELSEIF (ih(m04+i).EQ.'CH3 ') THEN !starts at i+1+1
799:                SEG_START(nseg+1) = i-1656:                SEG_START(nseg+1) = i-1
1025:       ELSE 882:       ELSE 
1026: 883: 
1027: C ******** backbone **********884: C ******** backbone **********
1028:          ! don't want to twist dihedrals in proline885:          ! don't want to twist dihedrals in proline
1029:          IF (RESLAB.EQ.'PRO'.OR.RESLAB.EQ.'NPRO'.OR.RESLAB.EQ.'CPRO') THEN886:          IF (RESLAB.EQ.'PRO'.OR.RESLAB.EQ.'NPRO'.OR.RESLAB.EQ.'CPRO') THEN
1030:             LSIDECHAIN = .FALSE.887:             LSIDECHAIN = .FALSE.
1031:             RETURN888:             RETURN
1032:          ENDIF889:          ENDIF
1033: 890: 
1034: ! ********** sidechain **************************************8891: ! ********** sidechain **************************************8
1035:          IF (RESLAB.EQ.'ARG'.OR.RESLAB.EQ.'NARG'.OR.RESLAB.EQ.'CARG')THEN892:          IF (RESLAB.EQ.'ASN'.OR.RESLAB.EQ.'NASN'.OR.RESLAB.EQ.'CASN') THEN
1036:             IF ((TYPEJ.EQ.'NE').AND.(TYPEK.EQ.'CZ')) LSIDECHAIN=.FALSE. 
1037:             IF ((TYPEJ.EQ.'CZ').AND.(TYPEK.EQ.'NH1')) LSIDECHAIN=.FALSE. 
1038:             IF ((TYPEJ.EQ.'CZ').AND.(TYPEK.EQ.'NH2')) LSIDECHAIN=.FALSE. 
1039:          ELSEIF (RESLAB.EQ.'ASN'.OR.RESLAB.EQ.'NASN'.OR.RESLAB.EQ.'CASN') THEN 
1040:             IF ((TYPEJ.EQ.'CG').AND.(TYPEK.EQ.'ND2')) LSIDECHAIN=.FALSE.893:             IF ((TYPEJ.EQ.'CG').AND.(TYPEK.EQ.'ND2')) LSIDECHAIN=.FALSE.
1041:             ! We do not twist the CG-ND2 bond since it is a peptide bond894:             ! We do not twist the CG-ND2 bond since it is a peptide bond
1042:          ELSEIF (reslab.EQ."HIE".OR.reslab.EQ.'NHIE'.OR.reslab.EQ.'CHIE') THEN895:          ELSEIF (reslab.EQ."HIE".OR.reslab.EQ.'NHIE'.OR.reslab.EQ.'CHIE') THEN
1043:             LSIDECHAIN = .FALSE.896:             LSIDECHAIN = .FALSE.
1044:             IF ((TYPEJ.EQ.'CA').AND.(TYPEK.EQ.'CB')) LSIDECHAIN=.TRUE.897:             IF ((TYPEJ.EQ.'CA').AND.(TYPEK.EQ.'CB')) LSIDECHAIN=.TRUE.
1045:             IF ((TYPEJ.EQ.'CB').AND.(TYPEK.EQ.'CG')) LSIDECHAIN=.TRUE.898:             IF ((TYPEJ.EQ.'CB').AND.(TYPEK.EQ.'CG')) LSIDECHAIN=.TRUE.
1046:             IF (.NOT.LSIDECHAIN)899:             IF (.NOT.LSIDECHAIN)
1047:      &       PRINT *, 'WARNING: His ring dihedrals will not be twisted.'900:      &       PRINT *, 'WARNING: His ring dihedrals will not be twisted.'
1048:          ELSEIF (RESLAB.EQ.'GLN'.OR.RESLAB.EQ.'NGLN'.OR.RESLAB.EQ.'CGLN') THEN901:          ELSEIF (RESLAB.EQ.'GLN'.OR.RESLAB.EQ.'NGLN'.OR.RESLAB.EQ.'CGLN') THEN
1049:             LSIDECHAIN = .FALSE.902:             LSIDECHAIN = .FALSE.
1089:      1      (TYPEK.EQ.'C').AND.(TYPEL.EQ.'N'))          .OR. 942:      1      (TYPEK.EQ.'C').AND.(TYPEL.EQ.'N'))          .OR. 
1090:      2       ((TYPEI.EQ.'N').AND.(TYPEJ.EQ.'CA').AND.943:      2       ((TYPEI.EQ.'N').AND.(TYPEJ.EQ.'CA').AND.
1091:      3       (TYPEK.EQ.'C').AND.(TYPEL.EQ.'NT'))         .OR.944:      3       (TYPEK.EQ.'C').AND.(TYPEL.EQ.'NT'))         .OR.
1092:      &       ((TYPEI.EQ.'N').AND.(TYPEJ.EQ.'CA')945:      &       ((TYPEI.EQ.'N').AND.(TYPEJ.EQ.'CA')
1093:      5        .AND.(TYPEK.EQ.'C').AND.(TYPEL.EQ.'OXT'))) THEN946:      5        .AND.(TYPEK.EQ.'C').AND.(TYPEL.EQ.'OXT'))) THEN
1094:              !PRINT *,'PSI',IICD 947:              !PRINT *,'PSI',IICD 
1095:          ENDIF948:          ENDIF
1096:       ENDIF949:       ENDIF
1097: 950: 
1098:       END SUBROUTINE ICTYPECHECKAM951:       END SUBROUTINE ICTYPECHECKAM
1099:  
1100:  
1101: ! ******************************************************************* 
1102:       SUBROUTINE GET_TWISTABLE(ST_PHI,FI_PHI) 
1103:       use modamber9 
1104:       IMPLICIT NONE 
1105: C This subroutine is important if perturbation is done in internal coordinates and  
1106: C if AMBPERTT is true 
1107: C calculates number of total twistable angles, NTW_ANGLES, sets elements which correspond 
1108: C to twistable coordinates in array TW_ANGLES to true and calculates P_DIFF, the  
1109: C probability of twisting according to difference between start and end structure 
1110:       DOUBLE PRECISION, INTENT(IN)::ST_PHI(nphia+nphih), FI_PHI(nphia+nphih) 
1111:       INTEGER :: A_SIDE, A_PHIPSI, A_OMEGA 
1112:       INTEGER :: i 
1113:       DOUBLE PRECISION :: DIFF_PHI(nphia+nphih) 
1114:       DOUBLE PRECISION :: MAX_DIFF 
1115:  
1116:       IF (AMBPERTT) THEN 
1117:          IF (.NOT. ALLOCATED(TW_DIFFP)) ALLOCATE(TW_DIFFP(lenic)) 
1118:          TW_DIFFP = -0.2D0 
1119:          NTW_ANGLES = 0 
1120:          A_SIDE = 1 
1121:          A_PHIPSI = 1 
1122:          A_OMEGA = 1 
1123:          DO i=1, lenic 
1124:             DIFF_PHI(i) = FI_PHI(i)-ST_PHI(i) 
1125:             IF (DIFF_PHI(i) .GT.180.0D0) DIFF_PHI(i) =  360.0D0 -DIFF_PHI(i) 
1126:             IF (DIFF_PHI(i) .LT.-180.0D0) DIFF_PHI(I) = 360.0D0 +DIFF_PHI(i) 
1127:          ENDDO 
1128:          MAX_DIFF = MAXVAL(DABS(DIFF_PHI)) 
1129:          DO i=1, lenic 
1130:             !take into account 2 pi circularity 
1131:             IF (PHIPSI(A_PHIPSI)==i) THEN 
1132:                IF (DABS(DIFF_PHI(i)).GT.PERTHRESH) THEN 
1133:                   TW_DIFFP(i) = (DABS(DIFF_PHI(i))/MAX_DIFF)/5.0D0 +0.20D0 
1134:                   !PRINT '(a15,i4,f7.3)',"twist phipsi",i, TW_DIFFP(i) 
1135:                ! to count how many twistable angles I have 
1136:                   NTW_ANGLES = NTW_ANGLES +1 
1137:                ENDIF 
1138:                A_PHIPSI = A_PHIPSI +1 
1139:             ELSEIF (TW_SIDECHAIN(A_SIDE)==i) THEN 
1140:                IF (DABS(DIFF_PHI(i)).GT.PERTHRESH) THEN 
1141:                   TW_DIFFP(i) =  (DABS(DIFF_PHI(i))/MAX_DIFF)/5.0D0 +0.20D0 
1142:                   !PRINT '(a15,i4,f7.3)',"twist side",i,TW_DIFFP(i)/5.0+0.2 
1143:                   NTW_ANGLES = NTW_ANGLES +1 
1144:                ENDIF 
1145:                A_SIDE = A_SIDE + 1 
1146:             ELSEIF (TOMEGAC) THEN 
1147:                IF (OMEGAC(A_OMEGA) ==i) THEN 
1148:                    IF (DABS(DIFF_PHI(i)).GT.PERTHRESH) THEN 
1149:                       TW_DIFFP(i) = DABS(DIFF_PHI(i))/MAX_DIFF 
1150:                       NTW_ANGLES = NTW_ANGLES +1 
1151:                    ENDIF 
1152:                    A_OMEGA = A_OMEGA +1 
1153:                ENDIF 
1154:             ENDIF 
1155:          ENDDO 
1156:       ENDIF 
1157:       PRINT*, "NTW_ANGLES", NTW_ANGLES 
1158:       END SUBROUTINE GET_TWISTABLE 


r14089/charmm_main.src 2017-01-21 10:32:28.079188235 +0000 r14088/charmm_main.src 2017-01-21 10:32:32.859188235 +0000
693:       ELSE IF (WRD.EQ.'HBON') THEN693:       ELSE IF (WRD.EQ.'HBON') THEN
694:         IF(IHBFRQ.EQ.0) IHBFRQ=999694:         IF(IHBFRQ.EQ.0) IHBFRQ=999
695:         CALL UPDATE(COMLYN,COMLEN,X,Y,Z,WMAIN,.TRUE.,695:         CALL UPDATE(COMLYN,COMLEN,X,Y,Z,WMAIN,.TRUE.,
696:      $        .FALSE.,.FALSE.,.TRUE.,.FALSE.,0,0,0,0,0,0,0)696:      $        .FALSE.,.FALSE.,.TRUE.,.FALSE.,0,0,0,0,0,0,0)
697:         IF(IHBFRQ.EQ.999) IHBFRQ=0697:         IF(IHBFRQ.EQ.999) IHBFRQ=0
698:       ELSE IF (WRD.EQ.'HBTR') THEN698:       ELSE IF (WRD.EQ.'HBTR') THEN
699:         CALL HBTRIM699:         CALL HBTRIM
700:       ELSE IF (WRD.EQ.'HBUI') THEN700:       ELSE IF (WRD.EQ.'HBUI') THEN
701:         CALL HBUILD(COMLYN,COMLEN)701:         CALL HBUILD(COMLYN,COMLEN)
702:       ELSE IF (WRD.EQ.'IC  ') THEN702:       ELSE IF (WRD.EQ.'IC  ') THEN
703:         PRINT*, "in chsetup l 691 IC (msb50)" 
704:         PRINT*, "msb50 complen", COMLEN 
705:         CALL INTCOR(COMLYN,COMLEN)703:         CALL INTCOR(COMLYN,COMLEN)
706:       ELSE IF (WRD.EQ.'INQU') THEN704:       ELSE IF (WRD.EQ.'INQU') THEN
707: C---- Procedure PROCESS-INQUIRE-FILE-COMMAND705: C---- Procedure PROCESS-INQUIRE-FILE-COMMAND
708: C706: C
709: C     inquire all open files707: C     inquire all open files
710:         IF(PRNLEV.GE.2) WRITE(OUTU,'(A)')708:         IF(PRNLEV.GE.2) WRITE(OUTU,'(A)')
711:      $        ' CHARMM: list of open files:'709:      $        ' CHARMM: list of open files:'
712:         DO I=1,99710:         DO I=1,99
713:           CALL VINQRE('UNIT',TNAME,TMAX,TLEN,QOPEN,QFORM,QWRITE,I)711:           CALL VINQRE('UNIT',TNAME,TMAX,TLEN,QOPEN,QFORM,QWRITE,I)
714:           IF (QOPEN) THEN712:           IF (QOPEN) THEN
831:       ELSE IF (WRD.EQ.'POLA') THEN829:       ELSE IF (WRD.EQ.'POLA') THEN
832:         CALL POLAR0830:         CALL POLAR0
833: ##ENDIF831: ##ENDIF
834: ##IF PATHINT832: ##IF PATHINT
835:       ELSE IF (WRD.EQ.'PINT') THEN833:       ELSE IF (WRD.EQ.'PINT') THEN
836:         CALL PINT(COMLYN,COMLEN)834:         CALL PINT(COMLYN,COMLEN)
837: ##ENDIF835: ##ENDIF
838:       ELSE IF (WRD.EQ.'PRES') THEN836:       ELSE IF (WRD.EQ.'PRES') THEN
839:         CALL GETPRS(COMLYN,COMLEN)837:         CALL GETPRS(COMLYN,COMLEN)
840:       ELSE IF (WRD.EQ.'PRIN') THEN838:       ELSE IF (WRD.EQ.'PRIN') THEN
841:         PRINT*, "in chsetup msb50 before mainio l827" 
842:         CALL MAINIO(WRD)839:         CALL MAINIO(WRD)
843:       ELSE IF (WRD.EQ.'PULL') THEN840:       ELSE IF (WRD.EQ.'PULL') THEN
844: C---- 23-Jul-96, LN841: C---- 23-Jul-96, LN
845:         CALL PULL(COMLYN,COMLEN)842:         CALL PULL(COMLYN,COMLEN)
846:       ELSE IF (WRD.EQ.'QUAN') THEN843:       ELSE IF (WRD.EQ.'QUAN') THEN
847:         CALL QMDEFN(COMLYN,COMLEN)844:         CALL QMDEFN(COMLYN,COMLEN)
848:       ELSE IF (WRD.EQ.'READ') THEN845:       ELSE IF (WRD.EQ.'READ') THEN
849:         CALL MAINIO(WRD)846:         CALL MAINIO(WRD)
850: ##IF GENETIC847: ##IF GENETIC
851: C---- 12-Jan-98, CLBIII848: C---- 12-Jan-98, CLBIII


r14089/Dijinit.f90 2017-01-21 10:32:32.359188235 +0000 r14088/Dijinit.f90 2017-01-21 10:32:36.979188235 +0000
177:             IF (TMPWEIGHT.GT.0.0D0) TMPWEIGHT=1.0D0/TMPWEIGHT177:             IF (TMPWEIGHT.GT.0.0D0) TMPWEIGHT=1.0D0/TMPWEIGHT
178:          ELSE178:          ELSE
179:             TMPWEIGHT=DUMMY**COSTFUNCTIONPOWER179:             TMPWEIGHT=DUMMY**COSTFUNCTIONPOWER
180:          ENDIF180:          ENDIF
181:       ENDIF181:       ENDIF
182:    ELSE182:    ELSE
183:       TMPWEIGHT=DUMMY183:       TMPWEIGHT=DUMMY
184:    ENDIF184:    ENDIF
185:    NSTEPS=NSTEPS+1185:    NSTEPS=NSTEPS+1
186:    186:    
187:    PRINT '(2(I8,G20.10),3G20.10)',J5,EMIN(J5),parent(J5),EMIN(PARENT(J5)),DUMMY/SCALEFAC,TMPWEIGHT,WEIGHT(J5)187:    PRINT '(2(I8,G20.10),3G20.10)',J5,EMIN(J5),parent(J5),EMIN(PARENT(J5)),DUMMY,TMPWEIGHT,WEIGHT(J5)
188:    THRESH=0.0D0188:    THRESH=0.0D0
189:    IF (BHINTERPT) THRESH=BHDISTTHRESH ! for bhinterp runs raise the threshold to BHDISTTHRESH189:    IF (BHINTERPT) THRESH=BHDISTTHRESH ! for bhinterp runs raise the threshold to BHDISTTHRESH
190:    IF (BISECTT) THRESH=BISECTMINDIST ! for bisect runs raise the threshold to BISECTMINDIST190:    IF (BISECTT) THRESH=BISECTMINDIST ! for bisect runs raise the threshold to BISECTMINDIST
191:    IF ((DUMMY/SCALEFAC.GT.THRESH).AND.(TMPWEIGHT.LT.HUGE(1.0D0)/10.0D0)) THEN191:    IF ((DUMMY/SCALEFAC.GT.THRESH).AND.(TMPWEIGHT.LT.HUGE(1.0D0)/10.0D0)) THEN
192:       NWORST=NWORST+1192:       NWORST=NWORST+1
193:       IF (NWORST.GT.10000) THEN193:       IF (NWORST.GT.10000) THEN
194:          PRINT '(A,I8)','ERROR in Dijinit, too many gaps, NWORST=',NWORST194:          PRINT '(A,I8)','ERROR in Dijinit, too many gaps, NWORST=',NWORST
195:          STOP195:          STOP
196:       ENDIF196:       ENDIF
197:       DMIN1(NWORST)=J5197:       DMIN1(NWORST)=J5


r14089/intcor2.src 2017-01-21 10:32:28.547188235 +0000 r14088/intcor2.src 2017-01-21 10:32:33.331188235 +0000
141:       INTEGER II1, JJ1, KK1, LL1141:       INTEGER II1, JJ1, KK1, LL1
142:       INTEGER QAT(5)142:       INTEGER QAT(5)
143:       INTEGER IX143:       INTEGER IX
144:       INTEGER NPERD144:       INTEGER NPERD
145:       REAL*8 EMIN,ANGMIN,EVAL,ANG145:       REAL*8 EMIN,ANGMIN,EVAL,ANG
146:       LOGICAL LDROP146:       LOGICAL LDROP
147:       INTEGER MARK2147:       INTEGER MARK2
148:       DATA MARK2/-99/148:       DATA MARK2/-99/
149: C149: C
150: C150: C
151:       PRINT*, "hello in Purgic" 
152:       IF(.NOT.LFILL) THEN151:       IF(.NOT.LFILL) THEN
153: C152: C
154:          PRINT*, "hello msb50 in intcor2" 
155:          IF(LCLEAR) THEN153:          IF(LCLEAR) THEN
156:             DO I=1,LENIC154:             DO I=1,LENIC
157:                IF(IAR(I).LE.0) IAR(I)=MARK2155:                IF(IAR(I).LE.0) IAR(I)=MARK2
158:                IF(JAR(I).LE.0) IAR(I)=MARK2156:                IF(JAR(I).LE.0) IAR(I)=MARK2
159:                IF(KAR(I).LE.0) IAR(I)=MARK2157:                IF(KAR(I).LE.0) IAR(I)=MARK2
160:                IF(LAR(I).LE.0) IAR(I)=MARK2158:                IF(LAR(I).LE.0) IAR(I)=MARK2
161:             ENDDO159:             ENDDO
162:          ENDIF160:          ENDIF
163: C161: C
164:          II=0162:          II=0
178:                   T1IC(II)=T1IC(I)176:                   T1IC(II)=T1IC(I)
179:                   T2IC(II)=T2IC(I)177:                   T2IC(II)=T2IC(I)
180:                   PIC(II)=PIC(I)178:                   PIC(II)=PIC(I)
181:                ENDIF179:                ENDIF
182:             ENDIF180:             ENDIF
183:          ENDDO181:          ENDDO
184:          LENIC=II182:          LENIC=II
185:       ENDIF183:       ENDIF
186: C184: C
187:       IF(LFILL) THEN185:       IF(LFILL) THEN
188:          PRINT*, "msb50 in l1081 intcor2.f" 
189:          PRINT*, "LENIC", LENIC 
190:          DO II=1,LENIC186:          DO II=1,LENIC
191:             IF(B1IC(II).LE.0.001 .OR.LALL) THEN187:             IF(B1IC(II).LE.0.001 .OR.LALL) THEN
192:                II1=KAR(II)188:                II1=KAR(II)
193:                JJ1=LAR(II)189:                JJ1=LAR(II)
194:                PRINT*, "msb50, B1IC, II1, JJ1", II1, JJ1 
195:                IF(II1.GT.0.AND.JJ1.GT.0) THEN190:                IF(II1.GT.0.AND.JJ1.GT.0) THEN
196:                   QAT(1)=II1191:                   QAT(1)=II1
197:                   QAT(2)=JJ1192:                   QAT(2)=JJ1
198:                   CALL CODES(QAT(3),0,0,0,193:                   CALL CODES(QAT(3),0,0,0,
199:      &                       NATOM,IMOVE,IAC,1,QAT(1),QAT(2),194:      &                       NATOM,IMOVE,IAC,1,QAT(1),QAT(2),
200:      &                       0,0,0,0,0,0,0,0,0,0,0,0,0,0,195:      &                       0,0,0,0,0,0,0,0,0,0,0,0,0,0,
201:      &                       .FALSE.,0,                   ! DRUDE196:      &                       .FALSE.,0,                   ! DRUDE
202:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP197:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP
203:      &                       .TRUE.,.TRUE.)198:      &                       .TRUE.,.TRUE.)
204:                   IX=QAT(3)199:                   IX=QAT(3)
205:                   PRINT*, "msb50 B1IC IX", IX 
206:                   IF(IX.GT.0) THEN200:                   IF(IX.GT.0) THEN
207:                      B1IC(II)=CBB(IX)201:                      B1IC(II)=CBB(IX)
208:                      PRINT*, "msb50 B1IC II",II, B1IC(II) 
209:                   ELSE202:                   ELSE
210:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(5A)')203:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(5A)')
211:      &               ' %BUILDC-WARNING: no bond parameters for types (',204:      &               ' %BUILDC-WARNING: no bond parameters for types (',
212:      &                ATC(IAC(II1)),',',ATC(IAC(JJ1)),')'205:      &                ATC(IAC(II1)),',',ATC(IAC(JJ1)),')'
213:                   ENDIF206:                   ENDIF
214:                ENDIF207:                ENDIF
215:             ENDIF208:             ENDIF
216: C        209: C        
217:             IF(B2IC(II).LE.0.001 .OR.LALL) THEN210:             IF(B2IC(II).LE.0.001 .OR.LALL) THEN
218:                II1=IAR(II)211:                II1=IAR(II)
221:                IF(II1.GT.0.AND.JJ1.GT.0) THEN214:                IF(II1.GT.0.AND.JJ1.GT.0) THEN
222:                   QAT(1)=II1215:                   QAT(1)=II1
223:                   QAT(2)=JJ1216:                   QAT(2)=JJ1
224:                   CALL CODES(QAT(3),0,0,0,217:                   CALL CODES(QAT(3),0,0,0,
225:      &                       NATOM,IMOVE,IAC,1,QAT(1),QAT(2),218:      &                       NATOM,IMOVE,IAC,1,QAT(1),QAT(2),
226:      &                       0,0,0,0,0,0,0,0,0,0,0,0,0,0,219:      &                       0,0,0,0,0,0,0,0,0,0,0,0,0,0,
227:      &                       .FALSE.,0,                   ! DRUDE220:      &                       .FALSE.,0,                   ! DRUDE
228:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP221:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP
229:      &                       .TRUE.,.TRUE.)222:      &                       .TRUE.,.TRUE.)
230:                   IX=QAT(3)223:                   IX=QAT(3)
231:                   PRINT*, "msb50, B2IC, IX", IX 
232:                   IF(IX.GT.0) THEN224:                   IF(IX.GT.0) THEN
233:                      B2IC(II)=CBB(IX)225:                      B2IC(II)=CBB(IX)
234:                      PRINT*, "msb50 B2IC II", II, B2IC(II) 
235:                   ELSE226:                   ELSE
236:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(5A)')227:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(5A)')
237:      &               ' %BUILDC-WARNING: no bond parameters for types (',228:      &               ' %BUILDC-WARNING: no bond parameters for types (',
238:      &                ATC(IAC(II1)),',',ATC(IAC(JJ1)),')'229:      &                ATC(IAC(II1)),',',ATC(IAC(JJ1)),')'
239:                   ENDIF230:                   ENDIF
240:                ENDIF231:                ENDIF
241:             ENDIF232:             ENDIF
242: C        233: C        
243:             IF(T1IC(II).LE.0.001 .OR.LALL) THEN234:             IF(T1IC(II).LE.0.001 .OR.LALL) THEN
244:                II1=LAR(II)235:                II1=LAR(II)
249:                   QAT(2)=JJ1240:                   QAT(2)=JJ1
250:                   QAT(3)=KK1241:                   QAT(3)=KK1
251:                   CALL CODES(0,QAT(4),0,0,242:                   CALL CODES(0,QAT(4),0,0,
252:      &                       NATOM,IMOVE,IAC,0,0,0,243:      &                       NATOM,IMOVE,IAC,0,0,0,
253:      &                       1,QAT(1),QAT(2),QAT(3),244:      &                       1,QAT(1),QAT(2),QAT(3),
254:      &                       0,0,0,0,0,0,0,0,0,0,245:      &                       0,0,0,0,0,0,0,0,0,0,
255:      &                       .FALSE.,0,                   ! DRUDE246:      &                       .FALSE.,0,                   ! DRUDE
256:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP247:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP
257:      &                       .TRUE.,.TRUE.)248:      &                       .TRUE.,.TRUE.)
258:                   IX=QAT(4)249:                   IX=QAT(4)
259:                   PRINT*, "msb50 T1IC,",II1,JJ1,KK1, IX 
260:                   IF(IX.GT.0) THEN250:                   IF(IX.GT.0) THEN
261:                      T1IC(II)=CTB(IX)*RADDEG251:                      T1IC(II)=CTB(IX)*RADDEG
262:                      PRINT*, "msb50 T1IC", T1IC(II) 
263:                   ELSE252:                   ELSE
264:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(7A)')253:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(7A)')
265:      &              ' %BUILDC-WARNING: no angle parameters for types (',254:      &              ' %BUILDC-WARNING: no angle parameters for types (',
266:      &             ATC(IAC(II1)),',',ATC(IAC(JJ1)),',',ATC(IAC(KK1)),')'255:      &             ATC(IAC(II1)),',',ATC(IAC(JJ1)),',',ATC(IAC(KK1)),')'
267:                   ENDIF256:                   ENDIF
268:                ENDIF257:                ENDIF
269:             ENDIF258:             ENDIF
270: C        259: C        
271:             IF(T2IC(II).LE.0.001 .OR.LALL) THEN260:             IF(T2IC(II).LE.0.001 .OR.LALL) THEN
272:                II1=IAR(II)261:                II1=IAR(II)
281:                   QAT(2)=JJ1270:                   QAT(2)=JJ1
282:                   QAT(3)=KK1271:                   QAT(3)=KK1
283:                   CALL CODES(0,QAT(4),0,0,272:                   CALL CODES(0,QAT(4),0,0,
284:      &                       NATOM,IMOVE,IAC,0,0,0,273:      &                       NATOM,IMOVE,IAC,0,0,0,
285:      &                       1,QAT(1),QAT(2),QAT(3),274:      &                       1,QAT(1),QAT(2),QAT(3),
286:      &                       0,0,0,0,0,0,0,0,0,0,275:      &                       0,0,0,0,0,0,0,0,0,0,
287:      &                       .FALSE.,0,                   ! DRUDE276:      &                       .FALSE.,0,                   ! DRUDE
288:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP277:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP
289:      &                       .TRUE.,.TRUE.)278:      &                       .TRUE.,.TRUE.)
290:                   IX=QAT(4)279:                   IX=QAT(4)
291:                   !PRINT*, "msb50 T2IC", II1,JJ1,KK1,IX 
292:                   IF(IX.GT.0) THEN280:                   IF(IX.GT.0) THEN
293:                      T2IC(II)=CTB(IX)*RADDEG281:                      T2IC(II)=CTB(IX)*RADDEG
294:                      PRINT*, "msb50 T2IC",T2IC(II)  
295:                   ELSE282:                   ELSE
296:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(7A)')283:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(7A)')
297:      &              ' %BUILDC-WARNING: no angle parameters for types (',284:      &              ' %BUILDC-WARNING: no angle parameters for types (',
298:      &             ATC(IAC(II1)),',',ATC(IAC(JJ1)),',',ATC(IAC(KK1)),')'285:      &             ATC(IAC(II1)),',',ATC(IAC(JJ1)),',',ATC(IAC(KK1)),')'
299:                   ENDIF286:                   ENDIF
300:                ENDIF287:                ENDIF
301:             ENDIF288:             ENDIF
302: C289: C
303: C           fill dihedrals (if explicitly requested)290: C           fill dihedrals (if explicitly requested)
304:             IF(LDIHE .AND. .NOT.TAR(II)) THEN291:             IF(LDIHE .AND. .NOT.TAR(II)) THEN
305:               PRINT*, "msb50, in dihedral run" 
306:               IF(PIC(II).LE.0.001 .OR.LALL) THEN292:               IF(PIC(II).LE.0.001 .OR.LALL) THEN
307:                 II1=IAR(II)293:                 II1=IAR(II)
308:                 JJ1=JAR(II)294:                 JJ1=JAR(II)
309:                 KK1=KAR(II)295:                 KK1=KAR(II)
310:                 LL1=LAR(II)296:                 LL1=LAR(II)
311:                 IF(II1.GT.0.AND.JJ1.GT.0.AND.KK1.GT.0.AND.LL1.GT.0) THEN297:                 IF(II1.GT.0.AND.JJ1.GT.0.AND.KK1.GT.0.AND.LL1.GT.0) THEN
312:                   QAT(1)=II1298:                   QAT(1)=II1
313:                   QAT(2)=JJ1299:                   QAT(2)=JJ1
314:                   QAT(3)=KK1300:                   QAT(3)=KK1
315:                   QAT(4)=LL1301:                   QAT(4)=LL1


r14089/intcor.src 2017-01-21 10:32:28.323188235 +0000 r14088/intcor.src 2017-01-21 10:32:33.087188235 +0000
 90: C 90: C
 91: C----------------------------------------------------------------------- 91: C-----------------------------------------------------------------------
 92: C 92: C
 93: C (flag for obvious coding errors - no reason to issue a message) 93: C (flag for obvious coding errors - no reason to issue a message)
 94:       IF(LENIC.GT.0  .AND.  LENIC.GT.INTLEN)  CALL DIEWRN(-4) 94:       IF(LENIC.GT.0  .AND.  LENIC.GT.INTLEN)  CALL DIEWRN(-4)
 95:       IF(LENICS.GT.0 .AND. LENICS.GT.INTLENS) CALL DIEWRN(-4) 95:       IF(LENICS.GT.0 .AND. LENICS.GT.INTLENS) CALL DIEWRN(-4)
 96: C 96: C
 97:       WRD=NEXTA4(COMLYN,COMLEN) 97:       WRD=NEXTA4(COMLYN,COMLEN)
 98: C 98: C
 99: C----------------------------------------------------------------------- 99: C-----------------------------------------------------------------------
100:       PRINT*, "hello in source" 
101:       IF(WRD.EQ.'READ') THEN100:       IF(WRD.EQ.'READ') THEN
102: C  READ-INTERNAL-COORDINATE-FILE101: C  READ-INTERNAL-COORDINATE-FILE
103:          IUNIC=GTRMI(COMLYN,COMLEN,'UNIT',-1)102:          IUNIC=GTRMI(COMLYN,COMLEN,'UNIT',-1)
104:          LAPPE=INDXA(COMLYN,COMLEN,'APPE').NE.0103:          LAPPE=INDXA(COMLYN,COMLEN,'APPE').NE.0
105:          ICARD=1104:          ICARD=1
106:          IF(INDXA(COMLYN,COMLEN,'FILE').NE.0) ICARD=0105:          IF(INDXA(COMLYN,COMLEN,'FILE').NE.0) ICARD=0
107:          IF(INDXA(COMLYN,COMLEN,'CARD').NE.0) ICARD=1106:          IF(INDXA(COMLYN,COMLEN,'CARD').NE.0) ICARD=1
108:          ISTART=1107:          ISTART=1
109:          IF(ICARD.EQ.1 .AND. IUNIC.LT.0) IUNIC=ISTRM108:          IF(ICARD.EQ.1 .AND. IUNIC.LT.0) IUNIC=ISTRM
110:          IF(INDXA(COMLYN,COMLEN,'SAVE').NE.0) THEN109:          IF(INDXA(COMLYN,COMLEN,'SAVE').NE.0) THEN
403:              I=LENIC+100  ! free excess space402:              I=LENIC+100  ! free excess space
404:              CALL REINTC(I,INTLEN,LENIC,BINTCR,LINTCR)403:              CALL REINTC(I,INTLEN,LENIC,BINTCR,LINTCR)
405:            ENDIF404:            ENDIF
406:          ENDIF405:          ENDIF
407: C-----------------------------------------------------------------------406: C-----------------------------------------------------------------------
408:       ELSE IF(WRD.EQ.'PARA') THEN407:       ELSE IF(WRD.EQ.'PARA') THEN
409: C  PROCESS-PARAMETER-COMMAND408: C  PROCESS-PARAMETER-COMMAND
410:          LAPPE=INDXA(COMLYN,COMLEN,'ALL').NE.0409:          LAPPE=INDXA(COMLYN,COMLEN,'ALL').NE.0
411:          LDIHE=INDXA(COMLYN,COMLEN,'DIHE').NE.0410:          LDIHE=INDXA(COMLYN,COMLEN,'DIHE').NE.0
412:          LIMPR=INDXA(COMLYN,COMLEN,'IMPR').NE.0411:          LIMPR=INDXA(COMLYN,COMLEN,'IMPR').NE.0
413:          PRINT*, "msb50 l413 intcor.src BINTCR" 
414:          PRINT*, BINTCR(INTB1),BINTCR(INTB2), BINTCR(INTT1) 
415:          PRINT*, BINTCR(INTIAR), BINTCR(INTJAR), BINTCR(INTKAR) 
416:          PRINT*, "msb50 l416 LENIC", LENIC 
417:          CALL PURGIC(LENIC,412:          CALL PURGIC(LENIC,
418:      &               HEAP(BINTCR(INTB1)),HEAP(BINTCR(INTB2)),413:      &               HEAP(BINTCR(INTB1)),HEAP(BINTCR(INTB2)),
419:      &               HEAP(BINTCR(INTT1)),HEAP(BINTCR(INTT2)),414:      &               HEAP(BINTCR(INTT1)),HEAP(BINTCR(INTT2)),
420:      &               HEAP(BINTCR(INTPIC)),HEAP(BINTCR(INTIAR)),415:      &               HEAP(BINTCR(INTPIC)),HEAP(BINTCR(INTIAR)),
421:      &               HEAP(BINTCR(INTJAR)),HEAP(BINTCR(INTKAR)),416:      &               HEAP(BINTCR(INTJAR)),HEAP(BINTCR(INTKAR)),
422:      &               HEAP(BINTCR(INTLAR)),HEAP(BINTCR(INTTAR)),417:      &               HEAP(BINTCR(INTLAR)),HEAP(BINTCR(INTTAR)),
423:      &               .TRUE.,LAPPE,LDIHE,LIMPR,.FALSE.,418:      &               .TRUE.,LAPPE,LDIHE,LIMPR,.FALSE.,
424:      &               NATOM,IMOVE,IAC)419:      &               NATOM,IMOVE,IAC)
425: C-----------------------------------------------------------------------420: C-----------------------------------------------------------------------
426:       ELSE IF(WRD.EQ.'DIFF') THEN421:       ELSE IF(WRD.EQ.'DIFF') THEN


r14089/keywords.f 2017-01-21 10:32:31.335188235 +0000 r14088/keywords.f 2017-01-21 10:32:35.935188235 +0000
 18: C 18: C
 19: C  All the keywords possible for the odata file are contained here in  19: C  All the keywords possible for the odata file are contained here in 
 20: C  alphabetical order. Initialisation statements procede the big IF block. 20: C  alphabetical order. Initialisation statements procede the big IF block.
 21: C 21: C
 22:       SUBROUTINE KEYWORDS(Q) 22:       SUBROUTINE KEYWORDS(Q)
 23:       USE COMMONS 23:       USE COMMONS
 24:       USE KEY 24:       USE KEY
 25:       USE MODMEC 25:       USE MODMEC
 26:       USE MODTWOEND 26:       USE MODTWOEND
 27:       USE MODAMBER 27:       USE MODAMBER
 28:       USE MODAMBER9, ONLY : COORDS1,ih,m04,prmtop,saltcon,igb,cut,rgbmax, 28:       USE MODAMBER9, ONLY : COORDS1,ih,m04,prmtop,saltcon,igb,cut,rgbmax,nocistransrna,atmass1,checkcistransalways, 
 29:      & nocistransrna,atmass1,checkcistransalways,AMBERICT,AMBSTEPT,AMBIT, 29:      &           ambpdb_unit, ambrst_unit, myunitnew, mdcrd_unit, mdinfo_unit
 30:      & AMBPERTT, PERTHRESH, AMBOLDPERTT,  
 31:      & ambpdb_unit, ambrst_unit, myunitnew, mdcrd_unit, mdinfo_unit 
 32:       USE MODNEB 30:       USE MODNEB
 33:       USE modmxatms   ! needed for CHARMM 31:       USE modmxatms   ! needed for CHARMM
 34:       USE modcharmm    32:       USE modcharmm   
 35:       USE MODUNRES 33:       USE MODUNRES
 36:       USE KEYNEB, NNNIMAGE=>NIMAGE 34:       USE KEYNEB, NNNIMAGE=>NIMAGE
 37:       USE KEYCONNECT 35:       USE KEYCONNECT
 38:       USE MODMEC 36:       USE MODMEC
 39:       USE MODGUESS 37:       USE MODGUESS
 40:       USE PORFUNCS 38:       USE PORFUNCS
 41:       USE msevb_common, ONLY: shellsToCount, maxHbondLength, minHbondAngle, OOclash_sq, printCoefficients 39:       USE msevb_common, ONLY: shellsToCount, maxHbondLength, minHbondAngle, OOclash_sq, printCoefficients
 57:  55: 
 58:       IMPLICIT NONE 56:       IMPLICIT NONE
 59:  57: 
 60:       INTEGER NDUM, NDIM1, NDIM2 58:       INTEGER NDUM, NDIM1, NDIM2
 61:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, IR, LAST, NTYPEA, J1, J2 59:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, IR, LAST, NTYPEA, J1, J2
 62:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR, 60:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR,
 63:      &                NERROR, IR, ECHO, LAST, CAT 61:      &                NERROR, IR, ECHO, LAST, CAT
 64:       DOUBLE PRECISION :: XX, EPSAB, EPSBB, SIGAB, SIGBB, Q(3*NATOMS), FDUM, RANDOM, DPRAND 62:       DOUBLE PRECISION :: XX, EPSAB, EPSBB, SIGAB, SIGBB, Q(3*NATOMS), FDUM, RANDOM, DPRAND
 65:       LOGICAL END, SKIPBL, CLEAR, ECHO, CAT, MINBFILE, CISTRANS 63:       LOGICAL END, SKIPBL, CLEAR, ECHO, CAT, MINBFILE, CISTRANS
 66:       CHARACTER WORD*25, WW*20 64:       CHARACTER WORD*25, WW*20
 67:       CHARACTER WORD2*25 
 68:       COMMON /BIN/ EPSAB, EPSBB, SIGAB, SIGBB, NTYPEA 65:       COMMON /BIN/ EPSAB, EPSBB, SIGAB, SIGBB, NTYPEA
 69:  66: 
 70:       INTEGER NATOM, DMODE 67:       INTEGER NATOM, DMODE
 71:       DOUBLE PRECISION CHX(MXATMS), CHY(MXATMS), CHZ(MXATMS), CHMASS(MXATMS) 68:       DOUBLE PRECISION CHX(MXATMS), CHY(MXATMS), CHZ(MXATMS), CHMASS(MXATMS)
 72:       DOUBLE PRECISION DPERT 69:       DOUBLE PRECISION DPERT
 73:       DOUBLE PRECISION CHPMIN, CHPMAX, CHNMIN, CHNMAX 70:       DOUBLE PRECISION CHPMIN, CHPMAX, CHNMIN, CHNMAX
 74:       INTEGER ISEED 71:       INTEGER ISEED
 75:  72: 
 76:       DOUBLE PRECISION UNRX(NATOMS), UNRY(NATOMS), UNRZ(NATOMS) ! UNRES 73:       DOUBLE PRECISION UNRX(NATOMS), UNRY(NATOMS), UNRZ(NATOMS) ! UNRES
 77:       DOUBLE PRECISION DUMMY1(NATOMS) 74:       DOUBLE PRECISION DUMMY1(NATOMS)
524:       BISECTDEBUG=.FALSE.521:       BISECTDEBUG=.FALSE.
525:       BISECTSTEPS=1522:       BISECTSTEPS=1
526:       BISECTMINDIST=1.0523:       BISECTMINDIST=1.0
527:       BISECTMAXATTEMPTS=5524:       BISECTMAXATTEMPTS=5
528:       CHRIGIDT=.FALSE.525:       CHRIGIDT=.FALSE.
529:       PTRANS=0.D0526:       PTRANS=0.D0
530:       TRANSMAX=0.D0527:       TRANSMAX=0.D0
531:       PROT=0.D0528:       PROT=0.D0
532:       ROTMAX=0.D0529:       ROTMAX=0.D0
533:       BBRSDMT=.FALSE.530:       BBRSDMT=.FALSE.
534:       AMBERICT=.FALSE.531:  
535:       AMBSTEPT=.FALSE. 
536:       AMBPERTT=.FALSE. 
537:       AMBOLDPERTT=.FALSE. 
538:       AMBIT = .FALSE.  
539:  
540:       DIJKSTRALOCAL=1.0D0532:       DIJKSTRALOCAL=1.0D0
541:       DNEBSWITCH=-1.0D0533:       DNEBSWITCH=-1.0D0
542:       PATHSDSTEPS=-1534:       PATHSDSTEPS=-1
543:       NUSEEV=-1535:       NUSEEV=-1
544: 536: 
545:       IF (FILTH2.EQ.0) THEN537:       IF (FILTH2.EQ.0) THEN
546:          OPEN (5,FILE='odata',STATUS='OLD')538:          OPEN (5,FILE='odata',STATUS='OLD')
547:       ELSE539:       ELSE
548:          WRITE(OTEMP,*) FILTH2540:          WRITE(OTEMP,*) FILTH2
549:          WRITE(OSTRING,'(A)') 'odata.' // TRIM(ADJUSTL(OTEMP))541:          WRITE(OSTRING,'(A)') 'odata.' // TRIM(ADJUSTL(OTEMP))
696: ! save atom names in array zsym688: ! save atom names in array zsym
697:         do i=1,natoms689:         do i=1,natoms
698:                 zsym(i) = ih(m04+i-1)690:                 zsym(i) = ih(m04+i-1)
699:         end do691:         end do
700:         RETURN692:         RETURN
701: ! initialise unit numbers693: ! initialise unit numbers
702:         ambpdb_unit=1110694:         ambpdb_unit=1110
703:         ambrst_unit=1111695:         ambrst_unit=1111
704:         mdinfo_unit=1112696:         mdinfo_unit=1112
705:         mdcrd_unit =1113697:         mdcrd_unit =1113
706:   
707:       ELSE IF (WORD.EQ.'AMBERIC') THEN 
708:         PRINT*, "amberic" 
709:         AMBERICT = .TRUE. 
710:         IF (NITEMS .GT. 1) THEN 
711:            CALL READA(WORD2) 
712:            IF (WORD2.EQ.'BACKBONE')  THEN 
713:               PRINT*, "backbone interpolated" 
714:               AMBIT = .TRUE. 
715:            ELSE  
716:               PRINT*, "keyword error in amberic" 
717:               RETURN 
718:            ENDIF 
719:         ENDIF 
720:   
721:       ELSE IF (WORD.eq.'AMBERSTEP') THEN 
722:         PRINT*, "amberstept" 
723:         AMBSTEPT = .TRUE.  
724:   
725:       ELSE IF (WORD.eq.'AMBPERTOLD') THEN 
726:         PRINT*, "original perturbation scheme" 
727:         AMBOLDPERTT = .TRUE. 
728:  
729:       ELSE IF (WORD.eq. 'AMBPERTONLY') THEN 
730:         AMBPERTT = .TRUE. 
731:         CALL READF(PERTHRESH)         
732:         PRINT*, "amber pertonly, perthresh", perthresh 
733:  
734:       ELSE IF (WORD.eq.'NAB') THEN698:       ELSE IF (WORD.eq.'NAB') THEN
735:         NABT=.TRUE.699:         NABT=.TRUE.
736:         DO i=1,3*NATOMS700:         DO i=1,3*NATOMS
737:                 Q(i) = COORDS1(i)701:                 Q(i) = COORDS1(i)
738:         END DO702:         END DO
739: ! save atom names in array zsym703: ! save atom names in array zsym
740:         do i=1,natoms704:         do i=1,natoms
741:                 zsym(i) = ih(m04+i-1)705:                 zsym(i) = ih(m04+i-1)
742:         end do706:         end do
743:         IF(.NOT.ALLOCATED(ATMASS)) ALLOCATE(ATMASS(NATOMS))707:         IF(.NOT.ALLOCATED(ATMASS)) ALLOCATE(ATMASS(NATOMS))
960: C924: C
961:       ELSE IF (WORD.EQ.'BSMIN') THEN925:       ELSE IF (WORD.EQ.'BSMIN') THEN
962:          BSMIN=.TRUE.926:          BSMIN=.TRUE.
963:          IF (NITEMS.GT.1) CALL READF(GMAX)927:          IF (NITEMS.GT.1) CALL READF(GMAX)
964:          IF (NITEMS.GT.2) CALL READF(EPS)928:          IF (NITEMS.GT.2) CALL READF(EPS)
965: 929: 
966:       ELSE IF (WORD.EQ.'BULK') THEN930:       ELSE IF (WORD.EQ.'BULK') THEN
967:          BULKT=.TRUE.931:          BULKT=.TRUE.
968: C932: C
969: C  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC933: C  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 934: C
970: C  CADPAC tells the program to read derivative information in 935: C  CADPAC tells the program to read derivative information in 
971: C         CADPAC format.                                        - default FALSE936: C         CADPAC format.                                        - default FALSE
972: C937: C
973:       ELSE IF (WORD.EQ.'CADPAC') THEN938:       ELSE IF (WORD.EQ.'CADPAC') THEN
974:          CADPAC=.TRUE.939:          CADPAC=.TRUE.
975:          CALL READA(SYS)940:          CALL READA(SYS)
976:          DO J1=1,80941:          DO J1=1,80
977:             IF (SYS(J1:J1).EQ.' ') THEN942:             IF (SYS(J1:J1).EQ.' ') THEN
978:                LSYS=J1-1943:                LSYS=J1-1
979:                GOTO 10944:                GOTO 10


r14089/modamber9.f90 2017-01-21 10:32:31.847188235 +0000 r14088/modamber9.f90 2017-01-21 10:32:36.251188235 +0000
297: integer, dimension(:), allocatable :: NSIDECHAIN297: integer, dimension(:), allocatable :: NSIDECHAIN
298: integer, dimension(:), allocatable :: NCHIRAL298: integer, dimension(:), allocatable :: NCHIRAL
299: integer, dimension(:), allocatable :: PHIPSI299: integer, dimension(:), allocatable :: PHIPSI
300: integer, dimension(:), allocatable :: OMEGAC300: integer, dimension(:), allocatable :: OMEGAC
301: integer, dimension(:), allocatable :: TW_SIDECHAIN301: integer, dimension(:), allocatable :: TW_SIDECHAIN
302: integer, dimension(:), allocatable :: CHIRAL302: integer, dimension(:), allocatable :: CHIRAL
303: logical, dimension(:), allocatable :: IS_SIDECHAIN303: logical, dimension(:), allocatable :: IS_SIDECHAIN
304: integer, dimension(:), allocatable :: NICTOT !no of residues per segment304: integer, dimension(:), allocatable :: NICTOT !no of residues per segment
305: integer:: nseg305: integer:: nseg
306: logical:: tomegac306: logical:: tomegac
307: logical:: AMBERICT, AMBSTEPT, AMBPERTT,AMBIT, AMBOLDPERTT 
308: double precision:: PERTHRESH !threshhold for dihedral perturbation 
309:                     !perturb if difference (start-end) > perthresh 
310: integer:: NTW_ANGLES !if ambpertt: how many angles can you really twist 
311: double precision, dimension(:), allocatable :: TW_DIFFP 
312: 307: 
313: 308: 
314: END MODULE modamber9 309: END MODULE modamber9 
315: 310: 


r14089/mylbfgs.f 2017-01-21 10:32:32.095188235 +0000 r14088/mylbfgs.f 2017-01-21 10:32:36.723188235 +0000
112:       IF (RESET) ITER=0112:       IF (RESET) ITER=0
113:       ITDONE=0113:       ITDONE=0
114:       IF (RESET.AND.PTEST) WRITE(*,'(A)') ' mylbfgs> Resetting LBFGS minimiser'114:       IF (RESET.AND.PTEST) WRITE(*,'(A)') ' mylbfgs> Resetting LBFGS minimiser'
115:       IF ((.NOT.RESET).AND.PTEST) WRITE(*,'(A)') ' mylbfgs> Not resetting LBFGS minimiser'115:       IF ((.NOT.RESET).AND.PTEST) WRITE(*,'(A)') ' mylbfgs> Not resetting LBFGS minimiser'
116: 1     FIXIMAGE=.FALSE.116: 1     FIXIMAGE=.FALSE.
117:       IF (PV.AND.(.NOT.PROJECT)) THEN117:       IF (PV.AND.(.NOT.PROJECT)) THEN
118:          IF (.NOT.KNOWE) THEN 118:          IF (.NOT.KNOWE) THEN 
119:             CALL POTENTIAL(X,ENERGY,GSAVE,.FALSE.,.FALSE.,RMS,.FALSE.,.FALSE.)119:             CALL POTENTIAL(X,ENERGY,GSAVE,.FALSE.,.FALSE.,RMS,.FALSE.,.FALSE.)
120: 120: 
121: ! Check for cold fusion121: ! Check for cold fusion
122:             PRINT*, "msb50 in mylbfgs, energy", energy 
123:             if (ENERGY.LT.coldFusionLimit) then122:             if (ENERGY.LT.coldFusionLimit) then
124:                ENERGY=1.0d60123:                ENERGY=1.0d60
125:                EREAL=1.0d60124:                EREAL=1.0d60
126:                RMS=1.0d1125:                RMS=1.0d1
127:                WRITE(*,'(A)') ' mylbfgs> Cold fusion diagnosed - step discarded'126:                WRITE(*,'(A)') ' mylbfgs> Cold fusion diagnosed - step discarded'
128:                RETURN127:                RETURN
129:             endif128:             endif
130:          ENDIF129:          ENDIF
131:          PVFLAG=.FALSE.130:          PVFLAG=.FALSE.
132:          CALL PVOPT(X,ENERGY,GSAVE)131:          CALL PVOPT(X,ENERGY,GSAVE)


r14089/OPTIM.tex 2017-01-21 10:32:28.919188235 +0000 r14088/OPTIM.tex 2017-01-21 10:32:33.603188235 +0000
153: starting coordinates.153: starting coordinates.
154: To turn on smooth cutoffs for the Generalised Born force fields, the keyword 154: To turn on smooth cutoffs for the Generalised Born force fields, the keyword 
155: {\it ifswitch=1} has to be used in the {\it \&cntrl} namelist block of {\it min.in}.155: {\it ifswitch=1} has to be used in the {\it \&cntrl} namelist block of {\it min.in}.
156: When using the {\it AMBER9} keyword, any calculated second derivatives will be 156: When using the {\it AMBER9} keyword, any calculated second derivatives will be 
157: numerical. If one wants analytical second derivatives, the {\it NAB} keyword 157: numerical. If one wants analytical second derivatives, the {\it NAB} keyword 
158: should be used instead, with the same syntax. The NAB interface does not 158: should be used instead, with the same syntax. The NAB interface does not 
159: currently support smooth cutoffs, so analytical second derivatives are only 159: currently support smooth cutoffs, so analytical second derivatives are only 
160: recommended for stationary points calculated with large cutoffs.160: recommended for stationary points calculated with large cutoffs.
161: Additional keywords for the AMBER 9 runs are {\it NAB} and {\it DUMPSTRUCTURES}.161: Additional keywords for the AMBER 9 runs are {\it NAB} and {\it DUMPSTRUCTURES}.
162: 162: 
163: \item{\it AMBERIC}: interpolation in bhinterp with internal coordinates. If BACKBONE specified additionally (AMBERIC BACKBONE) the protein backbone is also interpolated in internals. Does not work for polymers. 
164:  
165: \item{\it AMBERSTEP}: perturbation/twisting in bhinterp done with internal coordinates 
166:  
167: \item{\it AMBPERTONLY n}:only perturb angles for which the difference between the start and end configuration before the interpolation was greater than n (i.e. n= 30.0 - 30 degree). Probability of perturbation depends on this likelihood (and on position in chain) 
168:  
169: \item{\it AMBPERTOLD}: use original perturbation scheme. Now angles depend on position in chain - larger perturbation if at the end. 
170:  
171: \item {\it ANGLEAXIS}: specifies angle/axis coordinates for rigid body TIP potentials.163: \item {\it ANGLEAXIS}: specifies angle/axis coordinates for rigid body TIP potentials.
172: 164: 
173: \item {\it ARCTOL tol\/}: specifies the accuracy tolerance for computing the165: \item {\it ARCTOL tol\/}: specifies the accuracy tolerance for computing the
174:   inverse arc length used in the cubic spline interpolated string methods for166:   inverse arc length used in the cubic spline interpolated string methods for
175:   double-ended transition state searches. The default is $10^{-4}$.167:   double-ended transition state searches. The default is $10^{-4}$.
176: 168: 
177: \item {\it AXIS n}: specifies the highest symmetry axis to search for in169: \item {\it AXIS n}: specifies the highest symmetry axis to search for in
178: routine {\bf symmetry}; default is six.170: routine {\bf symmetry}; default is six.
179: 171: 
180: \item {\it BBCART\/}: use Cartesian coordinates for all backbone atoms and for172: \item {\it BBCART\/}: use Cartesian coordinates for all backbone atoms and for


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0