hdiff output

r15547/amber9dummy.f90 2017-01-21 10:33:11.490250000 +0000 r15546/amber9dummy.f90 2017-01-21 10:33:18.422250000 +0000
135: double precision :: THET2(1000), PHI(1000)135: double precision :: THET2(1000), PHI(1000)
136: logical :: DIHED_ONLY136: logical :: DIHED_ONLY
137: end subroutine CHGETICVAL137: end subroutine CHGETICVAL
138: 138: 
139: double precision function PHI_INTERP(ST_PHI,FI_PHI,SFRAC)139: double precision function PHI_INTERP(ST_PHI,FI_PHI,SFRAC)
140: double precision:: ST_PHI,FI_PHI,PHI,SFRAC140: double precision:: ST_PHI,FI_PHI,PHI,SFRAC
141: end function PHI_INTERP141: end function PHI_INTERP
142: 142: 
143: subroutine  TAKESTEPAMDIHED(COORDS, CSTART, CFINISH,SFRAC)143: subroutine  TAKESTEPAMDIHED(COORDS, CSTART, CFINISH,SFRAC)
144: double precision :: COORDS(3*1000), CSTART(3*1000), CFINISH(3*1000), SFRAC144: double precision :: COORDS(3*1000), CSTART(3*1000), CFINISH(3*1000), SFRAC
145: end subroutine TAKESTEPAMDIHED 
146: 145: 
147: SUBROUTINE AMBALIGNDIH(PHI,DIHC,RELPHI,PREVRELPHI,N)146: end subroutine TAKESTEPAMDIHED
148: double precision::PHI, RELPHI, PREVRELPHI, DRP, PHIORIG 
149: INTEGER:: DIHC, N 
150: end subroutine AMBALIGNDIH 
151:  
152: SUBROUTINE AMBGETOUTOFPLANE(I, J, K, L, PHI, dPdI, dPdJ, dPdK, dPdL, NOCOOR,NODERV) 
153: double precision:: I(0:2), J(0:2), K(0:2), L(0:2), PHI 
154: double precision:: dPdI(0:2), dPdJ(0:2), dPdK(0:2), dPdL(0:2) 
155: logical:: NOCOOR, NODERV 
156: end subroutine AMBGETOUTOFPLANE 
157:  
158: SUBROUTINE AMBGETANGLE(I, J, K, THETA, dTdI, dTdJ, dTdK, NOCOOR,NODERV) 
159: LOGICAL:: NOCOOR,NODERV 
160: DOUBLE PRECISION:: I(3), J(3), K(3), THETA 
161: DOUBLE PRECISION:: dTdI(3), dTdJ(3), dTdK(3) 
162: end subroutine AMBGETANGLE 
163:  
164: SUBROUTINE AMBGETTORSION(I, J, K, L, PHI, dPdI, dPdJ, dPdK, dPdL, NOCOOR, NODERV) 
165: LOGICAL:: NOCOOR, NODERV 
166: DOUBLE PRECISION:: I(3), J(3), K(3), L(3), PHI 
167: DOUBLE PRECISION:: dPdI(3), dPdJ(3), dPdK(3), dPdL(3) 
168: end subroutine AMBGETTORSION 
169:  
170: SUBROUTINE AMB_GETKDNAT(KD) 
171: INTEGER::KD 
172: end subroutine AMB_GETKDNAT 
173:  
174: subroutine AMBGETNNZNAT(NNZ) 
175: INTEGER::NNZ 
176: end subroutine AMBGETNNZNAT 
177:  
178: subroutine AMBTRANSFORM(XCART,GCART,XINT,GINT,NINTC,NCART,NNZ,NOCOOR,NODERV,KD,INTEPSILON) 
179: INTEGER:: NNZ,NINTC,NCART,KD 
180: DOUBLE PRECISION:: XCART(3*MXATMS),GCART(3*MXATMS),XINT(NINTC),GINT(NINTC) 
181: DOUBLE PRECISION:: INTEPSILON 
182: LOGICAL:: NOCOOR, NODERV 
183: end subroutine AMBTRANSFORM 
184:  
185: subroutine AMB_TRANSBACKDELTA(X,CARTX,COORDS,NINTC,NCART,NNZ,KD,FAILED,PTEST2,INTEPSILON) 
186: double precision INTEPSILON,CARTX(NCART),X(NINTC),COORDS(3*NATOM) 
187: integer NNZ,NINTC,NCART 
188: logical PTEST2,FAILEd 
189: end subroutine AMB_TRANSBACKDELTA 
190:  
191: subroutine check_valleu_chirality(COORDSB, COORDSA,DEBUG) 
192: double precision :: COORDSB(3*NATOMS), COORDSA(3*NATOMS) 
193: logical :: debug 
194: end subroutine check_valleu_chirality 
195:  
196: SUBROUTINE valleu_swap(COORDSB,COORDSA,i1,i4, DEBUG) 
197: double precision, intent(in):: COORDSB(3*NATOMS) 
198: double precision, intent(inout) :: COORDSA(3*NATOMS) 
199: integer, intent(in) ::i1,i4 
200: logical :: DEBUG 
201: END SUBROUTINE valleu_swap 
202: ======= 
203: 147: 
204: subroutine A9RESTOATOM(frozenres,frozen,nfreeze)148: subroutine A9RESTOATOM(frozenres,frozen,nfreeze)
205: integer i1, resnumber, currentresidue,nfreeze149: integer i1, resnumber, currentresidue,nfreeze
206: logical frozenres(natom), frozen(natom)150: logical frozenres(natom), frozen(natom)
207: 151: 
208: end subroutine A9RESTOATOM152: end subroutine A9RESTOATOM
209: 153: 
210: >>>>>>> .r15546 


r15547/amberinterface.f 2017-01-21 10:33:07.099188235 +0000 r15546/amberinterface.f 2017-01-21 10:33:14.610250000 +0000
1331:         angle=angle*180.0D0/pi1331:         angle=angle*180.0D0/pi
1332: !        write(*,*) 'calculated dihedral angle is', angle1332: !        write(*,*) 'calculated dihedral angle is', angle
1333:         IF (ABS(Q).LT.1.0D-10) PRINT '(A,G20.10)','calc_dihedral> WARNING, atoms nearly colinear, Q=',Q1333:         IF (ABS(Q).LT.1.0D-10) PRINT '(A,G20.10)','calc_dihedral> WARNING, atoms nearly colinear, Q=',Q
1334:       else1334:       else
1335:         write(*,*) 'ERROR - Q is 0! Are three atoms on the same line by any chance?'1335:         write(*,*) 'ERROR - Q is 0! Are three atoms on the same line by any chance?'
1336: !       STOP1336: !       STOP
1337:       end if1337:       end if
1338: 1338: 
1339: end subroutine calc_dihedral1339: end subroutine calc_dihedral
1340: 1340: 
 1341: 
1341:  subroutine cross2(u,v,p)1342:  subroutine cross2(u,v,p)
1342:  implicit none1343:  implicit none
1343:  double precision :: u(3),v(3),p(3)1344:  double precision :: u(3),v(3),p(3)
1344:  1345:  
1345:  p(1) = u(2)*v(3) - u(3)*v(2)1346:  p(1) = u(2)*v(3) - u(3)*v(2)
1346:  p(2) = u(3)*v(1) - u(1)*v(3)1347:  p(2) = u(3)*v(1) - u(1)*v(3)
1347:  p(3) = u(1)*v(2) - u(2)*v(1)1348:  p(3) = u(1)*v(2) - u(2)*v(1)
1348: 1349: 
1349:  end subroutine cross21350:  end subroutine cross2
1350: 1351: 
1556:         rij = SQRT(rij)1557:         rij = SQRT(rij)
1557:         IF (i+1==nbona) PRINT*, "now nbper"1558:         IF (i+1==nbona) PRINT*, "now nbper"
1558:         PRINT*, i3/3+1, j3/3+1, rij1559:         PRINT*, i3/3+1, j3/3+1, rij
1559:         bonds(nbonh +i+1) = rij1560:         bonds(nbonh +i+1) = rij
1560:       ENDDO1561:       ENDDO
1561:       END SUBROUTINE FINDALLCOORDS1562:       END SUBROUTINE FINDALLCOORDS
1562: 1563: 
1563: 1564: 
1564: ! ***********************************************************************1565: ! ***********************************************************************
1565: 1566: 
1566: subroutine check_valleu_chirality(COORDSB, COORDSA,DEBUG) 
1567: !msb50 - see whether order around prochiral VAL and LEU carbon is same in start and finish 
1568: !otherwise can't interpolate - and MINPERMDIST keeps swapping the groups if it thinks that's better 
1569:  use commons, only: NATOMS 
1570:  use modamber9 
1571:  implicit none 
1572:  double precision :: COORDSB(3*NATOMS), COORDSA(3*NATOMS) 
1573:  integer :: i,k, i1,i2,i3,i4 
1574:  double precision:: dihed1, dihed2 
1575:  double precision ::  TMP(3),TMP2(3),TMP3(3) 
1576:  logical :: debug   
1577:  
1578:  i1= -9999; i2= -9999; i3=-9999; i4=-9999 
1579:  dihed1 =0.0; dihed2 = 0.0 
1580:  do i = 1,nres 
1581:    IF ((ih(m02+i-1).EQ.'VAL').OR.(ih(m02+i-1).EQ.'NVAL').OR.(ih(m02+i-1) .EQ.'CVAL'))THEN 
1582:      DO k = ix(i02+i-1),ix(i02+i)-1 
1583:        IF (ih(m04+k-1).EQ.'CG1') i1 = k 
1584:        IF (ih(m04+k-1).EQ.'CA')  i2 = k 
1585:        IF (ih(m04+k-1).EQ.'CB')  i3 = k 
1586:        IF (ih(m04+k-1).EQ.'CG2') i4 = k 
1587:      ENDDO 
1588:      IF (i1.LT.0 .OR. i2.LT.0 .OR. i3.LT.0 .OR. i4.LT.0) THEN 
1589:          PRINT*, "check_valleu_chirality> VAL: ERROR" 
1590:          STOP 
1591:      ENDIF 
1592:      CALL AMBERDIHEDR(COORDSA, NATOMS,i1,i2,i3,i4, dihed1) 
1593:      CALL AMBERDIHEDR(COORDSB, NATOMS,i1,i2,i3,i4, dihed2) 
1594:      IF (dihed1*dihed2.LT.0) THEN !i.e. centres would have different chirality in start &end 
1595:         IF (DEBUG) PRINT*,"check_valleu_chirality: VAL",i,"chirality being changed" 
1596:         CALL valleu_swap(COORDSB,COORDSA,i1,i4,DEBUG)  
1597:      ENDIF 
1598:  
1599:    ELSEIF (ih(m02+i-1) .EQ.'LEU'.OR.ih(m02+i-1) .EQ.'NLEU'.OR.ih(m02+i-1) .EQ.'CLEU') THEN 
1600:      DO k = ix(i02+i-1),ix(i02+i)-1 
1601:        IF (ih(m04+k-1).EQ.'CD1') i1 = k 
1602:        IF (ih(m04+k-1).EQ.'CB')  i2 = k 
1603:        IF (ih(m04+k-1).EQ.'CG')  i3 = k 
1604:        IF (ih(m04+k-1).EQ.'CD2') i4 = k 
1605:      ENDDO 
1606:      IF (i1.LT.0 .OR. i2.LT.0 .OR. i3.LT.0 .OR. i4.LT.0) THEN 
1607:          PRINT*, "check_valleu_chirality> LEU: ERROR" 
1608:         STOP 
1609:      ENDIF 
1610:      CALL AMBERDIHEDR(COORDSB, NATOMS,i1,i2,i3,i4, dihed1) 
1611:      CALL AMBERDIHEDR(COORDSB, NATOMS,i1,i2,i3,i4, dihed2) 
1612:      IF (dihed1*dihed2.LT.0) THEN !i.e. centres would have different chirality in start &end 
1613:         IF (DEBUG) PRINT*,"check_valleu_chirality: LEU",i,"chirality being changed" 
1614:         CALL valleu_swap(COORDSB,COORDSA,i1,i4,DEBUG)  
1615:      ENDIF 
1616:  
1617:    ENDIF 
1618:  enddo 
1619: end subroutine check_valleu_chirality 
1620: ! ******************************* 
1621:  
1622: SUBROUTINE valleu_swap(COORDSB,COORDSA,i1,i4,debug) 
1623:  use commons, only: NATOMS 
1624:  use modamber9 
1625:  implicit none 
1626:  double precision, intent(in):: COORDSB(3*NATOMS) 
1627:  double precision, intent(inout) :: COORDSA(3*NATOMS) 
1628:  integer, intent(in) ::i1,i4  
1629:  integer ::ll 
1630:  logical :: debug 
1631:  double precision :: dihed1, dihed2, distance1, distance2 
1632:  double precision :: TMP(3), TMP2(3), TMP3(3), COORDS_copy(3*NATOMS) 
1633:   
1634:   COORDS_copy(:) = COORDSA(:) 
1635:   IF (i4-i1.EQ.4) THEN !all atoms 
1636:       COORDSA(3*(i1-1)+1:3*(i1+3-1)+3) = COORDS_copy(3*(i4-1)+1:3*(i4+3-1)+3) !swap carbon coords 
1637:       COORDSA(3*(i4-1)+1:3*(i4+3-1)+3) = COORDS_copy(3*(i1-1)+1:3*(i1+3-1)+3) 
1638:       COORDS_copy(:)= COORDSA(:) 
1639:       DO ll =i1,i4,i4-i1 !do for both 
1640:           CALL AMBERDIHEDR(COORDSA, NATOMS,ll+1,ll+2,ll,ll+3, dihed1) 
1641:           CALL AMBERDIHEDR(COORDSB, NATOMS,ll+1,ll+2,ll,ll+3, dihed2) 
1642:           IF (dihed1*dihed2.LT.0) THEN !wrong chirality - swap to H's  
1643:              COORDSA(3*(ll+1-1)+1:3*(ll+1-1)+3) = COORDS_copy(3*(ll+2-1)+1:3*(ll+2-1)+3) 
1644:              COORDSA(3*(ll+2-1)+1:3*(ll+2-1)+3) = COORDS_copy(3*(ll+1-1)+1:3*(ll+1-1)+3) 
1645:               distance1 = SQRT(DOT_PRODUCT(COORDSA, COORDSB)) 
1646:               TMP(:)= COORDSA(3*(ll)+1:3*(ll)+3) 
1647:               TMP2(:)= COORDSA(3*(ll+1)+1:3*(ll+1)+3) 
1648:               TMP3(:)= COORDSA(3*(ll+2)+1:3*(ll+2)+3) 
1649:               COORDSA(3*(ll-1):3*(ll-1)+3)=TMP2(:) 
1650:               COORDSA(3*(ll+1-1):3*(ll+1-1)+3) =TMP3(:) 
1651:               COORDSA(3*(ll+2-1):3*(ll+2-1)+3) =TMP(:) 
1652:               distance2 = SQRT(DOT_PRODUCT(COORDSA, COORDSB)) 
1653:               IF (distance1.LE.distance2) THEN !undo swap 
1654:                   COORDSA(:)= COORDS_copy(:) 
1655:               ELSE 
1656:                   distance1=distance2 
1657:                   IF (DEBUG) PRINT*, "permutation kept", ll+1, "->", ll+2, "etc" 
1658:               ENDIF 
1659:               COORDSA(3*(ll-1):3*(ll-1)+3)=TMP3(:) 
1660:               COORDSA(3*(ll+1-1):3*(ll+1-1)+3) =TMP(:) 
1661:               COORDSA(3*(ll+2-1):3*(ll+2-1)+3) =TMP2(:) 
1662:               distance2 = SQRT(DOT_PRODUCT(COORDSA, COORDSB)) 
1663:               IF (distance1.LE.distance2) THEN !undo swap 
1664:                   COORDSA(:)= COORDS_copy(:) 
1665:               ELSE 
1666:                   distance1=distance2 
1667:                   IF (DEBUG) PRINT*, "permutation kept", ll+1, "->", ll+3, "etc" 
1668:               ENDIF 
1669:          ENDIF 
1670:       ENDDO 
1671:   ELSEIF (i4-i1.EQ.1) THEN !united atoms 
1672:        COORDSA(3*(i1-1)+1:3*(i1-1)+3) = COORDS_copy(3*(i4-1)+1:3*(i4-1)+3) !swap carbon coords 
1673:        COORDSA(3*(i4-1)+1:3*(i4-1)+3) = COORDS_copy(3*(i1-1)+1:3*(i1-1)+3) 
1674:   ELSE 
1675:        PRINT*, "valleu_swap> Cannot determine force field. ERROR" 
1676:        STOP 
1677:   ENDIF 
1678: END SUBROUTINE valleu_swap 
1679:  
1680:  
1681:  
1682: !********************************* 
1683: 1567: 
1684:       SUBROUTINE TAKESTEPAMDIHED(Q, CSTART, CFINISH, SFRAC)1568:       SUBROUTINE TAKESTEPAMDIHED(Q, CSTART, CFINISH, SFRAC)
1685:       use modamber91569:       use modamber9
1686:       use commons, only : natoms1570:       use commons, only : natoms
1687:       IMPLICIT NONE1571:       IMPLICIT NONE
1688: ! similar to ICINTERP in twist.src1572: ! similar to ICINTERP in twist.src
1689: ! BOND1, THET1, PHI, THET2, BOND2 - IC arrays for interpolated structure1573: ! BOND1, THET1, PHI, THET2, BOND2 - IC arrays for interpolated structure
1690: ! GETICCOORDS generates the atoms defining the internal coordinates1574: ! GETICCOORDS generates the atoms defining the internal coordinates
1691: ! CHGETICVAL fills the internal coordinate table (similar to FILLICTABLE in Charmm)1575: ! CHGETICVAL fills the internal coordinate table (similar to FILLICTABLE in Charmm)
1692:       DOUBLE PRECISION, INTENT(IN) :: CSTART(3*natoms),CFINISH(3*natoms)1576:       DOUBLE PRECISION, INTENT(IN) :: CSTART(3*natoms),CFINISH(3*natoms)
1711:       DOUBLE PRECISION :: PI, RADDEG, DEGRAD1595:       DOUBLE PRECISION :: PI, RADDEG, DEGRAD
1712:       data pt999 /0.9990d0/1596:       data pt999 /0.9990d0/
1713:       PARAMETER(PI=3.141592653589793D0)1597:       PARAMETER(PI=3.141592653589793D0)
1714:       PARAMETER (RADDEG=180.0D0/PI)1598:       PARAMETER (RADDEG=180.0D0/PI)
1715:       PARAMETER (DEGRAD=PI/180.0D0)1599:       PARAMETER (DEGRAD=PI/180.0D0)
1716: !ask whether you have to consider the qpsander keyword 1600: !ask whether you have to consider the qpsander keyword 
1717: 1601: 
1718:       nterm_res =ix(i02+1)- ix(i02)1602:       nterm_res =ix(i02+1)- ix(i02)
1719: 1603: 
1720: 1604: 
1721:       PRINT*, "takestepamberdihed> ic"1605:       PRINT*, "IC"
1722:       IF (.NOT. ALLOCATED(IC_COORDS)) THEN1606:       IF (.NOT. ALLOCATED(IC_COORDS)) THEN
1723:          CALL GETICCOORDS()1607:          CALL GETICCOORDS()
1724:          PRINT*, "takestepamberdihed> lenic", lenic1608:          PRINT*, "lenic", lenic
1725:       ENDIF1609:       ENDIF
1726: 1610: 
1727:       IF (.NOT. ALLOCATED(IS_SIDECHAIN)) THEN1611:       IF (.NOT. ALLOCATED(IS_SIDECHAIN)) THEN
1728:           ALLOCATE(IS_SIDECHAIN(lenic))1612:           ALLOCATE(IS_SIDECHAIN(lenic))
1729:           DO i = 1, lenic1613:           DO i = 1, lenic
1730:               CALL CHECK_SIDECHAIN(IC_COORDS(i,1), IC_COORDS(i,2),      &1614:               CALL CHECK_SIDECHAIN(IC_COORDS(i,1), IC_COORDS(i,2),      &
1731:      &        IC_COORDS(i,3),IC_COORDS(i,4), i, IS_SIDECHAIN)1615:      &        IC_COORDS(i,3),IC_COORDS(i,4), i, IS_SIDECHAIN)
1732:           ENDDO1616:           ENDDO
1733:       ENDIF1617:       ENDIF
1734: 1618: 
2093:         rkj = vec_kj(1)*vec_kj(1)+vec_kj(2)*vec_kj(2)+vec_kj(3)*vec_kj(3)1977:         rkj = vec_kj(1)*vec_kj(1)+vec_kj(2)*vec_kj(2)+vec_kj(3)*vec_kj(3)
2094:         ct0 = vec_ij(1)*vec_kj(1)+vec_ij(2)*vec_kj(2)+vec_ij(3)*vec_kj(3)1978:         ct0 = vec_ij(1)*vec_kj(1)+vec_ij(2)*vec_kj(2)+vec_ij(3)*vec_kj(3)
2095:         ct0 = ct0/(SQRT(rij*rkj))1979:         ct0 = ct0/(SQRT(rij*rkj))
2096:         ct1 = max(-pt999,ct0) !still no idea why he's doing that1980:         ct1 = max(-pt999,ct0) !still no idea why he's doing that
2097:         ct2 = min(pt999,ct1)  !maybe making sure that ct2 != 11981:         ct2 = min(pt999,ct1)  !maybe making sure that ct2 != 1
2098:         ant = ACOS(ct2)1982:         ant = ACOS(ct2)
2099:         ant = ant*RADTOGRAD1983:         ant = ant*RADTOGRAD
2100:       END SUBROUTINE AMBERANGLE1984:       END SUBROUTINE AMBERANGLE
2101: 1985: 
2102: !**************1986: !**************
2103:  
2104:       SUBROUTINE GETICCOORDS()1987:       SUBROUTINE GETICCOORDS()
2105: 1988: 
2106:       use modamber91989:       use modamber9
2107:       use commons, only : NATOMS1990:       use commons, only : NATOMS
2108:       INTEGER :: i, j , a_index1991:       INTEGER :: i, j , a_index
2109:       INTEGER :: atomOffset,  nat_res, var1992:       INTEGER :: atomOffset,  nat_res, var
2110: !      INTEGER :: IC_COORDS(nres*95+5*nres,4)1993: !      INTEGER :: IC_COORDS(nres*95+5*nres,4)
2111:       INTEGER :: IC_START(23) !for each amino acid:where info in ICATOMS starts1994:       INTEGER :: IC_START(23) !for each amino acid:where info in ICATOMS starts
2112:       INTEGER :: IC_STARN(3), IC_STARC(3), IC_STARNPRO(3) !dummies for N/C terminals1995:       INTEGER :: IC_STARN(3), IC_STARC(3), IC_STARNPRO(3) !dummies for N/C terminals
2113:       CHARACTER(LEN=4) :: IC_ATOMS(1250)!(1121) !IMP if improper1996:       CHARACTER(LEN=4) :: IC_ATOMS(1250)!(1121) !IMP if improper
2363:      &                 "CE1", "CG", "CD1", "HD1", "IMP",                  &2246:      &                 "CE1", "CG", "CD1", "HD1", "IMP",                  &
2364:      &                 "CB", "CG", "CD2", "CE2", "NOT",                   &2247:      &                 "CB", "CG", "CD2", "CE2", "NOT",                   &
2365:      &                 "CE2", "CG", "CD2", "HD2", "IMP",                  &2248:      &                 "CE2", "CG", "CD2", "HD2", "IMP",                  &
2366:      &                 "CG", "CD1", "CE1", "CZ", "NOT",                   &2249:      &                 "CG", "CD1", "CE1", "CZ", "NOT",                   &
2367:      &                 "CZ", "CD1", "CE1", "HE1", "IMP",                  &2250:      &                 "CZ", "CD1", "CE1", "HE1", "IMP",                  &
2368:      &                 "CZ", "CD2", "CE2", "HE2", "IMP",                  &2251:      &                 "CZ", "CD2", "CE2", "HE2", "IMP",                  &
2369:      &                 "CE1", "CE2", "CZ", "OH", "IMP",                   &2252:      &                 "CE1", "CE2", "CZ", "OH", "IMP",                   &
2370:      &                 "CE1", "CZ", "OH", "HH", "NOT", &!end tyr           2253:      &                 "CE1", "CZ", "OH", "HH", "NOT", &!end tyr           
2371:      &                "N", "C", "CA", "CB", "IMP", &!val                   2254:      &                "N", "C", "CA", "CB", "IMP", &!val                   
2372:      &                "N", "C", "CA","HA", "IMP",                         &2255:      &                "N", "C", "CA","HA", "IMP",                         &
2373:      &                "N", "CA", "CB", "HB", "NOT",                      &2256:      &                "N", "CA", "CB", "CG1", "NOT",                      &
2374:      &                "HB", "CA", "CB", "CG1", "IMP",                    &2257:      &                "CG1", "CA", "CB", "CG2", "IMP",                    &
2375:      &                "HB", "CA", "CB", "CG2", "IMP",                    &2258:      &                "CG1", "CA", "CB", "HB", "IMP",                     &
2376: !     &                "N", "CA", "CB", "CG1", "NOT",                      & 
2377: !     &                "CG1", "CA", "CB", "CG2", "IMP",                    & 
2378: !     &                "CG1", "CA", "CB", "HB", "IMP",                     & 
2379:      &                "CA", "CB", "CG1", "HG11", "NOT",                   &2259:      &                "CA", "CB", "CG1", "HG11", "NOT",                   &
2380:      &                "HG11", "CB", "CG1", "HG12", "IMP",                 &2260:      &                "HG11", "CB", "CG1", "HG12", "IMP",                 &
2381:      &                "HG11", "CB", "CG1", "HG13", "IMP",                 &2261:      &                "HG11", "CB", "CG1", "HG13", "IMP",                 &
2382:      &                "CA", "CB", "CG2", "HG21", "NOT",                   &2262:      &                "CA", "CB", "CG2", "HG21", "NOT",                   &
2383:      &                "HG21", "CB", "CG2", "HG22", "IMP",                 &2263:      &                "HG21", "CB", "CG2", "HG22", "IMP",                 &
2384:      &                "HG21", "CB", "CG2", "HG23", "IMP"/ !end val        2264:      &                "HG21", "CB", "CG2", "HG23", "IMP"/ !end val        
2385: 2265: 
2386: 2266: 
2387: !      ix(i02...) number of first atoms in residues (where new resid starts)2267: !      ix(i02...) number of first atoms in residues (where new resid starts)
2388: !      ih(m02...) name of the residue2268: !      ih(m02...) name of the residue
2708:      &                   11 ,  6 ,  9 ,  12 ,0,                          &2588:      &                   11 ,  6 ,  9 ,  12 ,0,                          &
2709:      &                   4 ,  6 ,  10 ,  13 ,0,                          &2589:      &                   4 ,  6 ,  10 ,  13 ,0,                          &
2710:      &                   13 ,  6 ,  10 ,  14 ,0,                         &2590:      &                   13 ,  6 ,  10 ,  14 ,0,                         &
2711:      &                   6 ,  9 ,  11 ,  15 ,0,                          &2591:      &                   6 ,  9 ,  11 ,  15 ,0,                          &
2712:      &                   15 ,  9 ,  11 ,  16 ,0,                         &2592:      &                   15 ,  9 ,  11 ,  16 ,0,                         &
2713:      &                   15 ,  10 ,  13 ,  17 ,0,                        &2593:      &                   15 ,  10 ,  13 ,  17 ,0,                        &
2714:      &                   11 ,  13 ,  15 ,  18 ,0,                        &2594:      &                   11 ,  13 ,  15 ,  18 ,0,                        &
2715:      &                   11 ,  15 ,  18 ,  19 ,0,&!end tyr                &2595:      &                   11 ,  15 ,  18 ,  19 ,0,&!end tyr                &
2716:      &                   1 ,  2 ,  3 ,  4 ,0,&!val                        &2596:      &                   1 ,  2 ,  3 ,  4 ,0,&!val                        &
2717:      &                   1 ,  2 ,  3 ,  5 ,0,                            &2597:      &                   1 ,  2 ,  3 ,  5 ,0,                            &
2718:      &                   1,   3,   4 ,  6, 0,                            &2598:      &                   1 ,  3 ,  4 ,  6 ,0,                            &
2719:      &                   6 ,  3 ,  4 ,  7 ,0,                            &2599:      &                   6 ,  3 ,  4 ,  7 ,0,                            &
2720:      &                   6 ,  3 ,  4 ,  8 ,0,                            &2600:      &                   6 ,  3 ,  4 ,  8 ,0,                            &
2721:      &                   3 ,  4 ,  7 ,  9 ,0,                            &2601:      &                   3 ,  4 ,  6 ,  9 ,0,                            &
2722:      &                   9 ,  4 ,  7 ,  10 ,0,                           &2602:      &                   9 ,  4 ,  6 ,  10 ,0,                           &
2723:      &                   9 ,  4 ,  7 ,  11 ,0,                           &2603:      &                   9 ,  4 ,  6 ,  11 ,0,                           &
2724:      &                   3 ,  4 ,  8 ,  12 ,0,                           &2604:      &                   3 ,  4 ,  7 ,  12 ,0,                           &
2725:      &                   12 ,  4 ,  8 ,  13 ,0,                          &2605:      &                   12 ,  4 ,  7 ,  13 ,0,                          &
2726:      &                   12 ,  4 ,  8 ,  14 ,0/!end val2606:      &                   12 ,  4 ,  7 ,  14 ,0/!end val
2727: !     &                   3 ,  4 ,  6 ,  9 ,0,                            & 
2728: !     &                   9 ,  4 ,  6 ,  10 ,0,                           & 
2729: !     &                   9 ,  4 ,  6 ,  11 ,0,                           & 
2730: !     &                   3 ,  4 ,  7 ,  12 ,0,                           & 
2731: !     &                   12 ,  4 ,  7 ,  13 ,0,                          & 
2732: !     &                   12 ,  4 ,  7 ,  14 ,0/!end val 
2733: 2607: 
2734:         IF (IC_START(1) .EQ. 1) THEN !either NTERM or CTERM2608:         IF (IC_START(1) .EQ. 1) THEN !either NTERM or CTERM
2735:             var = IC_START(2)2609:             var = IC_START(2)
2736:             var2 = IC_START(3)-1 2610:             var2 = IC_START(3)-1 
2737:         ELSE2611:         ELSE
2738:            !could check len(IC_START ==21)      2612:            !could check len(IC_START ==21)      
2739:             var = IC_START(a_index)2613:             var = IC_START(a_index)
2740:             var2 = (IC_START(a_index+1)-var)/5-12614:             var2 = (IC_START(a_index+1)-var)/5-1
2741:         ENDIF2615:         ENDIF
2742: !      ENDIF2616: !      ENDIF


r15547/ambnatintern_extras.f90 2017-01-21 10:33:11.814250000 +0000 r15546/ambnatintern_extras.f90 2017-01-21 10:33:18.690250000 +0000
  1: ! **************************************************************************  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/ambnatintern_extras.f90' in revision 15546
  2:       
  3:       SUBROUTINE AMB_GETKDNAT(KD) 
  4: !get width of diagonal for GCS matrix, for natural internals 
  5: !for now, assume the dihedrals or bonds give the maximum width 
  6:       USE MODAMBER9 
  7:  
  8:  
  9:       INTEGER KD,ATOMWIDTH,DIFF,IPHI 
 10:       INTEGER I,J,K,L 
 11: ! 
 12:       ATOMWIDTH=0 
 13:       DIFF=0 
 14:       DO IPHI=1,NPHIA 
 15:         I=ix(i50+IPHI)/3 + 1 
 16:         J=ix(i52+IPHI)/3 + 1 
 17:         K=IABS(ix(i54+IPHI))/3 + 1 
 18:         L=IABS(ix(i56+IPHI))/3 + 1 
 19: !     print *,'pdihe I J K L',I,J,K,L 
 20:         DIFF=ABS(I-J) 
 21:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 22:         DIFF=ABS(I-K) 
 23:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 24:         DIFF=ABS(I-L) 
 25:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 26:         DIFF=ABS(J-K) 
 27:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 28:         DIFF=ABS(J-L) 
 29:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 30:         DIFF=ABS(K-L) 
 31:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 32:       ENDDO 
 33:       DO IPHI=1,NPHIH 
 34:         I=ix(i40+IPHI)/3 + 1 
 35:         J=ix(i42+IPHI)/3 + 1 
 36:         K=IABS(ix(i44+IPHI))/3 + 1 
 37:         L=IABS(ix(i46+IPHI))/3 + 1 
 38: !     print *,'pdihe I J K L',I,J,K,L 
 39:         DIFF=ABS(I-J) 
 40:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 41:         DIFF=ABS(I-K) 
 42:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 43:         DIFF=ABS(I-L) 
 44:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 45:         DIFF=ABS(J-K) 
 46:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 47:         DIFF=ABS(J-L) 
 48:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 49:         DIFF=ABS(K-L) 
 50:         IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 51:       ENDDO 
 52:  
 53:       IF(.NOT.ALLOCATED(IBIB)) THEN 
 54:          PRINT*, "ib going wrong!" 
 55:          ALLOCATE(IBIB(nbona+nbonh)) 
 56:          DO I =1,nbona 
 57:             IBIB(I) = ix(iiba+I-1)/3+1 
 58:          ENDDO 
 59:          DO I=1, nbonh 
 60:             IBIB(nbona+I) = ix(iibh+I-1)/3+1 
 61:          ENDDO 
 62:       ENDIF 
 63:  
 64:       IF(.NOT.ALLOCATED(JBJB)) THEN 
 65:           ALLOCATE(JBJB(nbona+nbonh)) 
 66:           DO I =1,nbona 
 67:              JBJB(I) = ix(ijba+I-1)/3+1 
 68:           ENDDO 
 69:           DO I=1, nbonh 
 70:              JBJB(nbona+I) = ix(ijbh+I-1)/3+1 
 71:           ENDDO 
 72:       ENDIF 
 73:  
 74:       DO IPHI = 1,nbona+nbonh 
 75:          I = IBIB(IPHI) 
 76:          J = JBJB(IPHI) 
 77:          DIFF=ABS(I-J) 
 78:          IF (DIFF.GT.ATOMWIDTH) ATOMWIDTH=DIFF 
 79:       ENDDO 
 80: ! 
 81:       KD=3*ATOMWIDTH+2 
 82: !     print *,'amb_getkdnat> ',KD 
 83: ! 
 84:       RETURN 
 85:       END 
 86:  
 87:  
 88: ! *************************************************************************** 
 89:  
 90:       SUBROUTINE AMBGETANGLE(I, J, K, THETA, dTdI, dTdJ, dTdK, NOCOOR,NODERV) 
 91: ! given the coordinates for three atoms I, J, K 
 92: ! find the angle with J as the central atom; return it in THETA 
 93: ! Also return all derivatives: dT/dIx, dT/dIy, dT/dIz in triplet dTdI 
 94: !same with dTdJ and dTdK 
 95:  
 96:       LOGICAL NOCOOR,NODERV 
 97:       DOUBLE PRECISION I(3), J(3), K(3), THETA 
 98:       DOUBLE PRECISION dTdI(3), dTdJ(3), dTdK(3) 
 99:       DOUBLE PRECISION DXI, DYI, DZI, DXJ, DYJ, DZJ, RI2, RJ2, RI, RJ, RIR, RJR 
100:       DOUBLE PRECISION DXIR, DYIR, DZIR, DXJR, DYJR, DZJR, CST, ISNT 
101:  
102:       DXI=I(1)-J(1) 
103:       DYI=I(2)-J(2) 
104:       DZI=I(3)-J(3) 
105:       DXJ=K(1)-J(1) 
106:       DYJ=K(2)-J(2) 
107:       DZJ=K(3)-J(3) 
108:       RI2=DXI*DXI+DYI*DYI+DZI*DZI 
109:       RJ2=DXJ*DXJ+DYJ*DYJ+DZJ*DZJ 
110:       RI=SQRT(RI2) 
111:       RJ=SQRT(RJ2) 
112:       RIR=1/RI 
113:       RJR=1/RJ 
114:       DXIR=DXI*RIR 
115:       DYIR=DYI*RIR 
116:       DZIR=DZI*RIR 
117:       DXJR=DXJ*RJR 
118:       DYJR=DYJ*RJR 
119:       DZJR=DZJ*RJR 
120:        
121:       CST=DXIR*DXJR+DYIR*DYJR+DZIR*DZJR 
122:       IF (.NOT.NOCOOR) THETA = ACOS(CST) 
123:  
124:       IF (.NOT.NODERV) THEN 
125:       ISNT=1/SQRT(1-CST*CST) 
126:       dTdI(1)=(DXIR*CST-DXJR)*RIR*ISNT 
127:       dTdI(2)=(DYIR*CST-DYJR)*RIR*ISNT 
128:       dTdI(3)=(DZIR*CST-DZJR)*RIR*ISNT 
129:       dTdJ(1)=(DXIR*(RI-RJ*CST)+DXJR*(RJ-RI*CST))*RIR*RJR*ISNT 
130:       dTdJ(2)=(DYIR*(RI-RJ*CST)+DYJR*(RJ-RI*CST))*RIR*RJR*ISNT 
131:       dTdJ(3)=(DZIR*(RI-RJ*CST)+DZJR*(RJ-RI*CST))*RIR*RJR*ISNT 
132:       dTdK(1)=(DXJR*CST-DXIR)*RJR*ISNT 
133:       dTdK(2)=(DYJR*CST-DYIR)*RJR*ISNT 
134:       dTdK(3)=(DZJR*CST-DZIR)*RJR*ISNT 
135:       ENDIF 
136:  
137:       RETURN 
138:       END SUBROUTINE 
139:  
140:       SUBROUTINE AMBGETOUTOFPLANE(I, J, K, L, PHI, dPdI, dPdJ, dPdK, dPdL, NOCOOR,NODERV) 
141: !EFK: 
142: !given the coordinates for four atoms I, J, K, L with I the central atom 
143: !find the angle between J and the plane defined by I,K,L; return it in PHI 
144: !Also return all derivatives: dP/dIx, dP/dIy, dP/dIz in triplet dPdI 
145: !same with dPdJ, dPdK, and dPdL 
146: !Calculations from Wilson, Delcius, and Cross (WDC) Ch.4 
147: ! 
148:       LOGICAL NOCOOR, NODERV 
149:       DOUBLE PRECISION I(0:2), J(0:2), K(0:2), L(0:2), PHI 
150:       DOUBLE PRECISION dPdI(0:2), dPdJ(0:2), dPdK(0:2), dPdL(0:2) 
151:       DOUBLE PRECISION JIX, JIY, JIZ, KIX, KIY, KIZ, LIX, LIY, LIZ 
152:       DOUBLE PRECISION RJI2, RKI2, RLI2, RJI, RKI, RLI, RJIR, RKIR, RLIR 
153:       DOUBLE PRECISION AX, AY, AZ, BX, BY, BZ, CX, CY, CZ 
154:       DOUBLE PRECISION RA2, RA, RB2, RB, RC2, RC, RCR 
155:       DOUBLE PRECISION SNP, CSP2, CSP, CSPR, TNP, DUMMY 
156:       INTEGER X 
157:  
158: !get vector differences 
159:       JIX=J(0)-I(0) ! r41 in WDC 
160:       JIY=J(1)-I(1) 
161:       JIZ=J(2)-I(2) 
162:       KIX=K(0)-I(0) ! r42 in WDC 
163:       KIY=K(1)-I(1) 
164:       KIZ=K(2)-I(2) 
165:       LIX=L(0)-I(0) ! r43 in WDC 
166:       LIY=L(1)-I(1) 
167:       LIZ=L(2)-I(2) 
168: !Get vector norms 
169:       RJI2 = JIX*JIX+JIY*JIY+JIZ*JIZ 
170:       RKI2 = KIX*KIX+KIY*KIY+KIZ*KIZ 
171:       RLI2 = LIX*LIX+LIY*LIY+LIZ*LIZ 
172:       RJI = SQRT(RJI2) 
173:       RKI = SQRT(RKI2) 
174:       RLI = SQRT(RLI2) 
175:       RJIR = 1/RJI 
176:       RKIR = 1/RKI 
177:       RLIR = 1/RLI 
178: !Get A = r41 x r42 
179:       AX = JIY*KIZ-KIY*JIZ 
180:       AY = -JIX*KIZ+KIX*JIZ 
181:       AZ = JIX*KIY-KIX*JIY 
182: !Get B = r43 x r41 
183:       BX = LIY*JIZ-JIY*LIZ 
184:       BY = -LIX*JIZ +JIX*LIZ 
185:       BZ = LIX*JIY-JIX*LIY 
186: !Get C = r42 x r43 
187:       CX = KIY*LIZ-LIY*KIZ 
188:       CY = -KIX*LIZ+LIX*KIZ 
189:       CZ = KIX*LIY-LIX*KIY 
190:       RC2 = CX*CX+CY*CY+CZ*CZ 
191:       RC = SQRT(RC2) 
192:       RCR = 1/RC 
193: !sin(phi1) = |r42 x r43| / (|r42| *|r43|) 
194: !DUMMY = (r42 x r43)*r41 
195:       DUMMY = CX*JIX+CY*JIY+CZ*JIZ 
196: !sin(PHI) = DUMMY / (sin(phi1) * |r42||r43||r41|) 
197:       SNP = DUMMY * RCR*RJIR  
198: !NOTE: POSSIBLE BUG; may run into problems if angle over 90 degrees 
199: !should fix this later 
200:       IF (.NOT.NOCOOR) PHI = ASIN(SNP) 
201:  
202:       IF (NODERV) RETURN 
203:  
204: !Calculate derivatives for B matrix 
205:       CSP2 = 1-SNP*SNP 
206:       CSP = SQRT(CSP2) !cos(PHI) 
207:       CSPR = 1/CSP !sec(PHI) 
208:       TNP = SNP*CSPR !tan(PHI) 
209: !POTENTIAL BUG: singularity at PHI->pi/2 
210: ! 
211:       dPdJ(0) = RJIR*(CX*CSPR*RCR - JIX*TNP*RJIR) 
212:       dPdJ(1) = RJIR*(CY*CSPR*RCR - JIY*TNP*RJIR) 
213:       dPdJ(2) = RJIR*(CZ*CSPR*RCR - JIZ*TNP*RJIR) 
214: ! 
215: !CONTINUE CHECKING WITH 2nd TERM HERE.... 
216:       DUMMY = LIX*KIX+LIY*KIY+LIZ*KIZ !r43*r42 
217:       dPdK(0) = BX*CSPR*RCR*RJIR - KIX*TNP*RLI2*RCR*RCR & 
218:      &     + LIX*TNP*DUMMY*RCR*RCR 
219:       dPdK(1) = BY*CSPR*RCR*RJIR - KIY*TNP*RLI2*RCR*RCR & 
220:      &     + LIY*TNP*DUMMY*RCR*RCR 
221:       dPdK(2) = BZ*CSPR*RCR*RJIR - KIZ*TNP*RLI2*RCR*RCR & 
222:      &     + LIZ*TNP*DUMMY*RCR*RCR 
223: !     
224:       dPdL(0) = AX*CSPR*RCR*RJIR - LIX*TNP*RKI2*RCR*RCR & 
225:      &     + KIX*TNP*DUMMY*RCR*RCR 
226:       dPdL(1) = AY*CSPR*RCR*RJIR - LIY*TNP*RKI2*RCR*RCR & 
227:      &     + KIY*TNP*DUMMY*RCR*RCR 
228:       dPdL(2) = AZ*CSPR*RCR*RJIR - LIZ*TNP*RKI2*RCR*RCR & 
229:      &     + KIZ*TNP*DUMMY*RCR*RCR 
230: ! 
231:       DO X = 0,2 
232:          dPdI(X) = -dPdJ(X)-dPdK(X)-dPdL(X) 
233:       ENDDO 
234: ! 
235:       RETURN 
236:       END 
237:  
238:       SUBROUTINE AMBGETTORSION(I, J, K, L, PHI, dPdI, dPdJ, dPdK, dPdL, NOCOOR, NODERV) 
239: !given the coordinates for four atoms I, J, K, L, bound in order,  
240: !find the dihedral torsion angle; return it in PHI 
241: !Also return all derivatives: dP/dIx, dP/dIy, dP/dIz in triplet dPdI 
242: !same with dPdJ, dPdK, and dPdL 
243:  
244:       LOGICAL NOCOOR, NODERV 
245:       DOUBLE PRECISION I(3), J(3), K(3), L(3), PHI 
246:       DOUBLE PRECISION dPdI(3), dPdJ(3), dPdK(3), dPdL(3) 
247:       DOUBLE PRECISION FX, FY, FZ, GX, GY, GZ, HX, HY, HZ 
248:       DOUBLE PRECISION AX, AY, AZ, BX, BY, Bz 
249:       DOUBLE PRECISION RF, RG, RH, RF2, RG2, RH2, RFR, RGR, RHR 
250:       DOUBLE PRECISION CSTTWO, SNTTWO2, CSTTHREE, SNTTHREE2, SNTTWO2R, SNTTHREE2R 
251:       DOUBLE PRECISION RA2, RB2, RA2R, RB2R, RABR, CP 
252:       DOUBLE PRECISION MYTX, MYTY, MYTZ, MYSCALAR 
253:       DOUBLE PRECISION DUMMY, DUMMY2 
254:  
255:       FX=I(1)-J(1) 
256:       FY=I(2)-J(2) 
257:       FZ=I(3)-J(3) 
258:       GX=J(1)-K(1) 
259:       GY=J(2)-K(2) 
260:       GZ=J(3)-K(3) 
261:       HX=L(1)-K(1) 
262:       HY=L(2)-K(2) 
263:       HZ=L(3)-K(3) 
264: !A=F^G, B=H^G 
265:       AX=FY*GZ-FZ*GY 
266:       AY=FZ*GX-FX*GZ 
267:       AZ=FX*GY-FY*GX 
268:       BX=HY*GZ-HZ*GY 
269:       BY=HZ*GX-HX*GZ 
270:       BZ=HX*GY-HY*GX 
271: !RG=|G|, RGR=1/|G| 
272:       RG2=GX*GX+GY*GY+GZ*GZ 
273:       RG=SQRT(RG2) 
274:       RGR=1/RG 
275: !dae for use in evaluating B-matrix 
276:       RF2=FX*FX+FY*FY+FZ*FZ 
277:       RF=SQRT(RF2) 
278:       RFR=1/RF 
279:       RH2=HX*HX+HY*HY+HZ*HZ 
280:       RH=SQRT(RH2) 
281:       RHR=1/RH 
282: !    jmc        CSTTWO=(FX*GX+FY*GY+FZ*GZ)*RFR*RGR 
283:       CSTTWO=-(FX*GX+FY*GY+FZ*GZ)*RFR*RGR 
284:       SNTTWO2=1-CSTTWO*CSTTWO 
285:       SNTTWO2R=1/SNTTWO2 
286:       CSTTHREE=(HX*GX+HY*GY+HZ*GZ)*RHR*RGR 
287:       SNTTHREE2=1-CSTTHREE*CSTTHREE 
288:       SNTTHREE2R=1/SNTTHREE2       
289: ! 
290:       IF (.NOT.NOCOOR) THEN 
291:          RA2=AX*AX+AY*AY+AZ*AZ 
292:          RB2=BX*BX+BY*BY+BZ*BZ 
293:          RA2R=1/RA2 
294:          RB2R=1/RB2 
295:          RABR=SQRT(RA2R*RB2R) 
296: !    CP=cos(phi) 
297:          CP=(AX*BX+AY*BY+AZ*BZ)*RABR 
298:          IF (CP.GT.1.0D0) THEN 
299:             PHI=0.0D0 
300:          ELSEIF (CP.LT.-1.0D0) THEN 
301:             PHI=PI 
302:          ELSE 
303:             PHI=ACOS(CP) 
304:          ENDIF 
305:  
306: !    jmc Test to determine whether the dihedral angle should be > 0 or < 0, from unres (intcor.f) 
307: !    > 0 if cw rotation from 2-> 1 to 3-> 4. 
308:          MYTX=AY*BZ-BY*AZ 
309:          MYTY=-AX*BZ+BX*AZ 
310:          MYTZ=AX*BY-BX*AY 
311:          MYSCALAR=MYTX*GX+MYTY*GY+MYTZ*GZ 
312:          IF (MYSCALAR.GT.0.0D0) PHI = -PHI 
313:       ENDIF 
314:  
315:       IF (NODERV) RETURN 
316:  
317: !dae calculate B-matrix elements 
318:       DUMMY=RFR*RFR*RGR*SNTTWO2R 
319:       dPdI = (/-AX*DUMMY, -AY*DUMMY, -AZ*DUMMY/) 
320:  
321:       DUMMY=RFR*RFR*RGR*RGR*SNTTWO2R*(RG-RF*CSTTWO) 
322:       DUMMY2=RHR*RGR*RGR*SNTTHREE2R*CSTTHREE 
323:       dPdJ = (/AX*DUMMY-BX*DUMMY2, AY*DUMMY-BY*DUMMY2, AZ*DUMMY-BZ*DUMMY2/) 
324:  
325:       DUMMY=RHR*RHR*RGR*RGR*SNTTHREE2R*(RG-RH*CSTTHREE) 
326:       DUMMY2=RFR*RGR*RGR*SNTTWO2R*CSTTWO 
327:       dPdK = (/-BX*DUMMY+AX*DUMMY2, -BY*DUMMY+AY*DUMMY2,-BZ*DUMMY+AZ*DUMMY2/) 
328:  
329:       DUMMY=RHR*RHR*RGR*SNTTHREE2R 
330:       dPdL = (/BX*DUMMY,BY*DUMMY,BZ*DUMMY/) 
331:  
332:       RETURN 
333:       END SUBROUTINE 
334:  
335:  
336: !************************************************************************************************ 
337:  
338:       SUBROUTINE AMBALIGNDIH(PHI,DIHC,RELPHI,PREVRELPHI,N) 
339: !aligning dihedral angles properly 
340: !if aligndir is set, and it's not the first dihedral in the group (N.NE.1) then  
341: !adjust angle to minimise change in relative dihedral to the first in the group 
342: !otherwise adjust angle to minimise change from previous value of same angle 
343:  
344:       USE INTCOMMONS, ONLY : PREVDIH, ALIGNDIR 
345:  
346:       IMPLICIT NONE 
347:  
348:       DOUBLE PRECISION PHI, RELPHI, PREVRELPHI, DRP, PHIORIG 
349:       INTEGER DIHC, N 
350:       DOUBLE PRECISION PI, TWOPI 
351:       PARAMETER(PI=3.141592653589793D0) 
352:  
353:       TWOPI = 2.0D0*PI 
354:       PHIORIG = PHI 
355:  
356:       IF (ALIGNDIR.AND.N.NE.1) THEN 
357:          DRP = RELPHI-PREVRELPHI 
358:          IF (DRP.GT.PI) THEN 
359:             PHI = PHI - TWOPI 
360:          ELSE IF (DRP.LT.-PI) THEN 
361:             PHI = PHI + TWOPI 
362:          ENDIF          
363:       ELSE 
364:          IF (PHI-PREVDIH(DIHC).GT.PI) THEN 
365:             PHI = PHI - TWOPI 
366:          ELSE IF (PHI-PREVDIH(DIHC).LT.-PI) THEN 
367:             PHI = PHI + TWOPI 
368:          ENDIF 
369:       ENDIF 
370:   
371:       PREVDIH(DIHC) = PHI 
372:  
373:       RETURN 
374:       END SUBROUTINE 
375:  
376:  
377: !********************************************************************* 
378:       SUBROUTINE AMBGETNNZNAT(NNZ) 
379: !    get number of nonzero elements in B matrix, for natural internals       
380:       USE COMMONS, ONLY : NATOMS 
381:       USE INTCOMMONS, ONLY : CARTATMSTART, ATOMxPOSinC, CENTERS, NBDS,& 
382:      &     NIMP, NFRG, NCRT, NCNT, LINDIH, RINGS, NRNG, NLDH 
383:  
384:       IMPLICIT NONE 
385:  
386:       INTEGER NNZ, I, NANGC, CNTTYPE, NATMSC 
387:  
388:       NNZ = 6*NBDS+12*NIMP+18*NFRG + 3*NCRT 
389:       DO I=1,NCNT 
390:          CNTTYPE = CENTERS(I,0) 
391:          IF (CNTTYPE.LT.3) THEN  
392:             NNZ = NNZ + ATOMxPOSinC(CNTTYPE,5,5) + 3 
393:          ELSE IF (CNTTYPE.EQ.3) THEN 
394:             NNZ = NNZ + ATOMxPOSinC(CNTTYPE,5,4) + 3 
395:          ELSE IF (CNTTYPE.EQ.4.OR.CNTTYPE.EQ.5) THEN 
396:             NNZ = NNZ + ATOMxPOSinC(CNTTYPE,2,4) + 3 
397:          ELSE IF (CNTTYPE.EQ.6) THEN 
398:             NNZ = NNZ + ATOMxPOSinC(CNTTYPE,1,3) + 3 
399:          ELSE IF (CNTTYPE.EQ.7) THEN 
400:             NNZ = NNZ + ATOMxPOSinC(CNTTYPE,1,4) + 3 
401:          ELSE IF (CNTTYPE.EQ.8) THEN 
402:             NNZ = NNZ +ATOMxPOSinC(CNTTYPE,4,5) +3 
403: !    ELSE IF (CNTTYPE.EQ.8) THEN 
404: !    NNZ = NNZ + ATOMxPOSinC(3,4,4) + 3 
405: !    ELSE IF (CNTTYPE.EQ.9) THEN 
406: !    NNZ = NNZ + ATOMxPOSinC(2,4,5) + 3 
407:          ENDIF 
408:       ENDDO 
409:       DO I=1,NRNG 
410:          IF(RINGS(I,6).EQ.-1) THEN 
411:             NNZ = NNZ + 3*5*4 
412:          ELSE 
413:             NNZ = NNZ + 3*6*6 
414:          ENDIF 
415:       ENDDO 
416:       DO I=1,NLDH 
417:          NNZ = NNZ + 6 + 3*(LINDIH(I,-1)+LINDIH(I,0)) 
418:       ENDDO 
419: !     print*, 'getnnznat>> nonzeros found is:', NNZ 
420:       RETURN 
421:       END 
422:  
423:  
424: !*********************************************************************** 
425: ! 
426: !subroutine to do internal coordinate manipulations using B-matrix 
427: !calls GETBEE to calcualte internal coordinates and B-matrix 
428: !then transforms cartesian gradients to internal gradients 
429: !Nemeth et al. JCP 113,5598(2000) 
430: ! 
431:  
432:       SUBROUTINE AMBTRANSFORM(XCART,GCART,XINT,GINT,NINTC,NCART,NNZ,NOCOOR,NODERV,KD,INTEPSILON) 
433:       USE modmxatms, ONLY: MXATMS 
434:   
435:        USE SPFUNCTS, ONLY : DUMPCOORDS 
436:  
437:        IMPLICIT NONE 
438: !SAT: two implicit none - line below commented! 
439: !#INCLUDE '~/charmmcode/fcm/dimens.fcm' 
440: !#INCLUDE '~/charmmcode/fcm/psf.fcm' 
441: !NNZ is number of non-zero elements in B-matrix 
442: !NINTC is number of internal coordinates 
443: !NCART is number of cartesians 
444: ! 
445:        INTEGER NNZ,NINTC,NCART,KD,K 
446: ! 
447:        DOUBLE PRECISION XCART(3*MXATMS),GCART(3*MXATMS),XINT(NINTC),GINT(NINTC) 
448:        DOUBLE PRECISION VAL(NNZ) 
449:        INTEGER INDX(NNZ),JNDX(NNZ) 
450:        INTEGER TRANSA,M,N,DESCRA(9),I1,J1,K1,LDGCS,LDA,LDB,INFO,NRHS,LWORK,ATOMWIDTH 
451:        DOUBLE PRECISION ALPHA,BETA,WORK(NINTC),GCS(KD+1,NCART),GCS2(2*KD+1,NCART) 
452:        DOUBLE PRECISION INTEPSILON,GX(NCART),GY(NINTC),SUM,SUMGY2,MEANSUMGY2,GCSREAL 
453:        CHARACTER*1 UPLO 
454:        LOGICAL NOCOOR,LCONVERGE, NODERV 
455:        DOUBLE PRECISION :: GCARTORIG(NCART), GNORM, GCARTINPUT(NCART) 
456:        INTEGER :: NITS 
457:  
458: !--------------------------------- 
459:  
460: !get B-matrix and internal coordiantes 
461: ! 
462:  
463:        CALL AMB_GETBEENAT(XCART,XINT,VAL,INDX,JNDX,NINTC,NCART, NNZ,GCS2,NOCOOR,NODERV,KD) 
464:        IF (NODERV) RETURN 
465:  
466:        GINT(1:NINTC) = 0.0D0 
467: !     print *,'tf here1' 
468: ! 
469: !VAL is column of non-zero BEE elements indexed by INDX and JNDX 
470: !GCS2 = BEEtBEE as calculated in GETBEE, stored in sparse format 
471: !need to convert double sparse format of GCS2 to lapack upperdiagonal only GCS 
472: !KD is number of superdiagonals of GC 
473: !see dpbtrf.f http://www.netlib.org/lapack/double/dpbtrf.f 
474: !      
475: !move to array of dimension KD+1 
476:        DO J1=1,KD 
477:           DO I1=1,J1 
478:              GCS(KD+1+I1-J1,J1)=GCS2(KD+1+I1-J1,J1) 
479:           ENDDO 
480:        ENDDO 
481:        DO J1=KD+1,NCART 
482:           DO I1=J1-KD,J1 
483:              GCS(KD+1+I1-J1,J1)=GCS2(KD+1+I1-J1,J1) 
484:           ENDDO 
485:        ENDDO 
486:  
487: ! 
488: !XINT now contains internal coordinates 
489: !now need to transform gradients 
490: ! 
491: !form GC = GC + eI where e is small 
492: ! 
493:        DO I1=1,NCART 
494:           GCS(KD+1,I1)=GCS(KD+1,I1)+INTEPSILON 
495:        ENDDO 
496:  
497: ! 
498: !Cholesky decompose GCS 
499: !use lapack routine DPBTRF 
500: !LDGCS is leading dimension of GCS 
501: ! 
502:           UPLO='U' 
503:           N=NCART 
504:           LDGCS=KD+1 
505:  
506:           CALL DPBTRF(UPLO,N,KD,GCS,LDGCS,INFO) 
507:  
508:           IF (INFO.NE.0) PRINT*,'WARNING - after DPBTRF INFO=',INFO 
509: !     print *,'tf here2' 
510: ! 
511: !GCS is now the upper triangular matrix after decomposition 
512: !SOLVE G*GCART2=GCART FOR GCART2, WHERE G IS ORIGINAL GC(=BEETBEE)+INTEPSILON*I 
513: !done using lapack routine DPBTRS with input of the upper matrix from 
514: !cholesky decomposition 
515: ! 
516:           NRHS=1 
517:           LDB=NCART 
518:           CALL DPBTRS(UPLO,N,KD,NRHS,GCS,LDGCS,GCART,LDB,INFO) 
519:           IF (INFO.NE.0) PRINT*,'WARNING - after DPBTRS INFO=',INFO 
520:  
521: !    print *,'tf here3' 
522: ! 
523: !GCART is now GCART2 
524: !GINT = BEE*GCART as first approximation 
525: !use DCOOMM to do this sparse multiplication 
526: !sparse matrix multiplication routine from BLAS 
527: !http://math.nist.gov/~KRemington/fspblas/ 
528: ! 
529:        TRANSA=0 
530:        M=NINTC 
531:        N=1 
532:        K=NCART 
533:        ALPHA=1.0D0 
534:        BETA=0.0D0 
535:        DO I1=1,9 
536:          DESCRA(I1)=0 
537:        ENDDO 
538:          DESCRA(4)=1 
539:        LWORK=NINTC 
540: !     print *,'tf here3.5' 
541:        CALL DCOOMM(TRANSA,M,N,K,ALPHA,DESCRA,VAL,INDX,JNDX,NNZ,GCART,NCART,BETA,GINT,NINTC,WORK,LWORK) 
542:        IF (INFO.NE.0) PRINT*,'WARNING - after DCOOMM INFO=',INFO 
543:  
544: !      IF(ITERATEG) PREITGRAD(:) = GINT(:) 
545: !     print *,'tf here4' 
546: ! 
547: !can skip below iteration and use this first approximation 
548: ! 
549: !print out GINT 
550: !     DO I1=1,NINTC 
551: !      print *,'GINT',GINT(I1) 
552: !     ENDDO 
553: !C 
554: !now iterate GINT 
555: !GINT(K+1)=GINT(K)+At(GRAD-BEEt*GINT(K)) 
556: !multiplication At*GX is the same as solving G*GCART2=GX 
557: !C 
558: !     print *,'transform here4' 
559: !      LCONVERGE = .FALSE. 
560: !      NITS = 0 
561: !      IF (ITERATEG) THEN 
562: !         DO WHILE (.NOT.LCONVERGE) 
563: !    GX=BEEt*GINT 
564: !GX=GCART-GX 
565: !            TRANSA=1 
566: !            M=NINTC; N = 1; K = NCART 
567: !            ALPHA=1.0D0 
568: !            BETA=0.0D0 
569: !            DO I1=1,9 
570: !               DESCRA(I1)=0 
571: !            ENDDO 
572: !            DESCRA(4)=1 
573: !            LWORK=NINTC 
574: !            GX(:) = 0.0D0 
575: !            CALL DCOOMM(TRANSA,M,N,K,ALPHA,DESCRA,VAL,INDX,JNDX,NNZ, 
576: !    $            GINT,NINTC,BETA,GX,NCART,WORK,LWORK) 
577: !            DO I1=1,NCART 
578: !               GX(I1)=GCARTORIG(I1)-GX(I1) 
579: !            ENDDO 
580: !C!             print*, 'TESTXC: ', SQRT(DOT_PRODUCT(GX,GX)/NCART) 
581: !C!GY=At(GRAD-BEEt*GINT(K))=At*GX 
582: !            NRHS=1 
583: !            LDB=NCART 
584: !            UPLO='U' 
585: !            N=NCART 
586: !            LDGCS=KD+1 
587: !            CALL DPBTRS(UPLO,N,KD,NRHS,GCS,LDGCS,GX,LDB,INFO) 
588: !            IF (INFO.NE.0) PRINT*,'WARNING - after DPBTRS 2 INFO=',INFO 
589: !C!       CALL DPOTRS(UPLO,N,NRHS,GC,LDA,GX,LDB,INFO) 
590: !            TRANSA=0 
591: !            M=NINTC; N = 1; K=NCART 
592: !            GY(:) = 0.0D0 
593: !            CALL DCOOMM(TRANSA,M,N,K,ALPHA,DESCRA,VAL,INDX,JNDX,NNZ,GX, 
594: !    $            NCART,BETA,GY,NINTC,WORK,LWORK)  
595:  
596: !             print*, 'TESTXI: ', SQRT(DOT_PRODUCT(GY,GY)/NINTC) 
597: !update GINT  
598: !            SUMGY2=0.0 
599: !            DO I1=1,NINTC 
600: !               GINT(I1)=GINT(I1)+GY(I1) 
601: !               SUMGY2=SUMGY2+GY(I1)*GY(I1) 
602: !            ENDDO 
603: !            MEANSUMGY2=SQRT(SUMGY2/NINTC) 
604: !            IF (MEANSUMGY2.LT.ITGRADCUTOFF) THEN 
605: !               LCONVERGE=.TRUE. 
606: !            ENDIF 
607: !             
608: !            NITS = NITS+1 
609: !            IF (NITS.GT.1000) THEN 
610: !               print*, 'ERROR: Grad iteration failed to converge!' 
611: !               CALL DUMPCOORDS(XCART,'badcoords.xyz', .FALSE.) 
612: !               CALL DUMPCOORDS(GCARTORIG, 'badgrad.xyz', .FALSE.) 
613: !               STOP 
614: !            ENDIF 
615: !      print *,'transform MEANSUMGY2',MEANSUMGY2 
616: !         ENDDO 
617: !      ENDIF 
618: !C 
619: !move final internal gradients to GRAD 
620: !C 
621: !     DO I1=1,NINTC 
622: !       GRAD(I1)=GINT(I1) 
623: !     ENDDO 
624: !C 
625: !zero GCS2 for next call to getbee 
626: !     DO J1=1,KD 
627: !       DO I1=1,J1 
628: !         GCS2(KD+1+I1-J1,J1)=0.0 
629: !         GCS2(KD+1+J1-I1,I1)=0.0 
630: !       ENDDO 
631: !     ENDDO 
632: !     DO J1=KD+1,NCART 
633: !       DO I1=J1-KD,J1 
634: !         GCS2(KD+1+I1-J1,J1)=0.0 
635: !         GCS2(KD+1+J1-I1,I1)=0.0 
636: !       ENDDO 
637: !     ENDDO 
638: !C 
639:        RETURN 
640:        END 
641: ! 
642:  
643:  
644: ! ************************************************************************* 
645:  
646: !      
647: !     subroutine to do internal coordinate manipulations using B-matrix 
648: !     calls GETBEE to calcualte B-matrix 
649: !     then transforms a step in internals to a step in cartesians  
650: !     Nemeth et al. JCP 113,5598(2000) 
651: !      
652:       SUBROUTINE AMB_TRANSBACKDELTA(X,CARTX,COORDS,NINTC,NCART,NNZ,KD,FAILED,PTEST2,INTEPSILON) 
653:       USE COMMONS, only: NATOMS 
654:       USE MODAMBER9, only: nbona, nbonh, ntheth, ntheta 
655:       USE INTCOMMONS, only : INTNEWT, NATINT, BACKTCUTOFF, BACKTPRINTERR 
656:       USE SPFUNCTS, only : DUMPCOORDS 
657:  
658:       IMPLICIT NONE 
659: !      
660:  
661: !     REAL*8 PREVDIH 
662: !     COMMON /DIHINFO/ PREVDIH(MAXP) 
663: !      
664: !     NNZ is number of non-zero elements in B-matrix 
665: !     NINTC is number of internal coordinates 
666: !     NCART is number of cartesians 
667: !    E  
668:       INTEGER NNZ,NINTC,NCART 
669: !      
670: !     COORDS is cartesian coordinates 
671: !     OLDQ is old internal coordinates 
672: !     X is eigenvector in internals, CARTX is ev in cartesians which we need to calculate 
673: !     need to zero CARTX 
674: !      
675:       DOUBLE PRECISION COORDS(3*NATOMS),DUMCOORDS(3*NATOMS) 
676:       DOUBLE PRECISION VAL(NNZ),WORK(NINTC) 
677:       INTEGER INDX(NNZ),JNDX(NNZ),NITS 
678:       DOUBLE PRECISION ALPHA,BETA 
679:       INTEGER TRANSA,M,N,DESCRA(9),I1,LDGCS,LDB,INFO,NRHS,LWORK,ATOMWIDTH,KD 
680:       DOUBLE PRECISION INTEPSILON,CARTX(NCART),X(NINTC),NEWX(NINTC),XNEW2(NINTC),DUMX(NINTC),SUMDC2,DUMMYQ(NINTC) 
681:       DOUBLE PRECISION MEANDC2,GCS2(2*KD+1,NCART),GCS(KD+1,NCART) 
682:       DOUBLE PRECISION CARTDOTX,CARTDOTY,CARTDOTZ 
683:       CHARACTER*1 UPLO 
684:       LOGICAL NOCOOR,LCONVERGE,FAILED,PTEST, NODERV 
685:       INTEGER J1, J2, K 
686: !    EFK 
687:       INTEGER DUMINDX(NNZ), DUMJNDX(NNZ) 
688:       DOUBLE PRECISION DUMVAL(NNZ), DUMGCS2(2*KD+1,NCART), DUMXINT1(NINTC), DUMXINT2(NINTC) 
689:       DOUBLE PRECISION QORIG(NINTC), TESTQ(NINTC), QTARGET(NINTC) 
690:       LOGICAL PTEST2 
691:  
692:  
693: !      INTEGER CENTERS, RINGS, FRINGS, LINDIH, IMPDIH, NCNT, NRNG, NFRG, 
694: !     $     NLDH, NIMP, NBDS, ATOMxPOSINC 
695: !      REAL*8 COEFF 
696: !      COMMON /NATINTERN/ COEFF(30,6,6), 
697: !     $     CENTERS(MAXA, 0:5), RINGS(MAXA, 6), 
698: !     $     FRINGS(MAXA, 6), LINDIH(MAXA, -1:8), IMPDIH(MAXA, 4), NCNT, 
699: !     $     NRNG, NFRG, NLDH, NIMP, NBDS, ATOMxPOSINC(10,5,5) 
700:        
701:       CHARACTER*30 FNAME 
702: !msb50 
703:       INTEGER NBOND, NTHETS 
704:       DOUBLE PRECISION PI, TWOPI 
705:       PARAMETER(PI=3.141592653589793D0) 
706:  
707:  
708:       TWOPI = 2.0D0*PI 
709:       NBOND = nbona+nbonh  
710:       NTHETS = ntheta +ntheth 
711:  
712: !      PTEST=.TRUE. 
713:       PTEST = .FALSE. 
714:  
715: !     make dummy array NEWX that can change (error in estimate of X) 
716: !     and DUMX which is cumulative estimate of X 
717: !     and dummy array DUMCOORDS that can change 
718: !      
719:       DO I1=1,NINTC 
720:          NEWX(I1)=X(I1) 
721:          DUMX(I1)=0.0D0 
722:       ENDDO 
723:       DO I1=1,NCART 
724:          CARTX(I1)=0.0D0 
725:          DUMCOORDS(I1)=COORDS(I1) 
726:       ENDDO 
727: !     
728: !    start iterative loop to do conversion of internal coordinates to cartesians 
729: !    have X, need to satisfy CARTX=A*X and X=B*CARTX 
730: !     
731:       LCONVERGE=.FALSE. 
732:       FAILED=.TRUE. 
733:       NITS=0 
734: !     
735: !    Attempt to prevent NaN on Solaris 
736: !     
737:       DUMMYQ(1:NINTC) = 0.0D0 
738:       XNEW2(1:NINTC)=0.0D0 
739:  
740:       DO WHILE (.NOT.LCONVERGE) 
741: !        IF (PTEST) print '(A,I6,G20.10)', 
742: !    $        'transbackdelta>> internals error', NITS, 
743: !    $        SQRT(DOT_PRODUCT(NEWX(1:NINTC),NEWX(1:NINTC))) 
744: !     
745: !    get B-matrix without changing cartesians to internals 
746: !    added DUMMYQ argument 13/2/04 DJW 
747: !     
748:          IF (PTEST) THEN 
749:             PRINT*,'DUMCOORDS BEFORE GETBEE, INTEPSILON=',INTEPSILON 
750:             WRITE(*,'(6G18.10)') (DUMCOORDS(J1),J1=1,NCART) 
751:             PRINT*,'DUMMYQ before GETBEE' 
752:             WRITE(*,'(6G18.10)') (DUMMYQ(J1),J1=1,NINTC) 
753:          ENDIF 
754:  
755:          IF (NITS.EQ.0) THEN ! store original and target internals 
756:             NOCOOR = .FALSE.; NODERV = .FALSE. 
757:             CALL AMB_GETBEENAT(DUMCOORDS,QORIG,VAL,INDX,JNDX,NINTC,NCART,NNZ,GCS2,NOCOOR,NODERV,KD) 
758:             QTARGET(:) = QORIG(:) + X(:) 
759:          ELSE 
760:             NOCOOR = .TRUE.; NODERV = .FALSE. 
761:             CALL AMB_GETBEENAT(DUMCOORDS,DUMMYQ,VAL,INDX,JNDX,NINTC,NCART,NNZ,GCS2,NOCOOR,NODERV,KD) 
762:          ENDIF 
763:  
764:         IF (PTEST) THEN 
765:            PRINT*,'GCS2 AFTER GETBEE:',INTEPSILON 
766:            WRITE(*,'(6G18.10)') ((GCS2(J1,J2),J1=1,2*KD+1),J2=1,NCART) 
767:         ENDIF 
768:  
769: !     
770: !    want to find CARTX=A*X 
771: !    A = G^-1*BEEt 
772: !    same as solving G*CARTX=BEEt*X 
773: !    G=GC+INTEPSILON*I 
774: !    GC=BEEtBEE (as in TRANSFORM), stored as sparse GCS 
775: !    KD is number of superdiagonals of GC 
776: !    see dpbtrf.f http://www.netlib.org/lapack/double/dpbtrf.f 
777: !     
778: !    KD=3*ATOMWIDTH+2        
779: !    need to move GCS2 into GCS 
780:  
781:          GCS(:,:) = 0.0D0 
782:          DO J1=1,KD 
783:             DO I1=1,J1 
784:                GCS(KD+1+I1-J1,J1)=GCS2(KD+1+I1-J1,J1) 
785:             ENDDO 
786:          ENDDO 
787:          DO J1=KD+1,NCART 
788:             DO I1=J1-KD,J1 
789:                GCS(KD+1+I1-J1,J1)=GCS2(KD+1+I1-J1,J1) 
790:                !print*, GCS(KD+1+I1-J1,J1) 
791:             ENDDO 
792:          ENDDO 
793:             
794: !     
795: !    form G = GC + eI where e is small 
796: !     
797:          DO I1=1,NCART 
798:             GCS(KD+1,I1)=GCS(KD+1,I1)+INTEPSILON 
799:          ENDDO 
800: !     
801: !    NEWX is error in internals at end of last cycle 
802: !    multiply NEWX by BEEt to give CARTX 
803: !    Uses sparse matrix multiplication routine from BLAS 
804: !    http://math.nist.gov/~KRemington/fspblas/ 
805: !     
806:          TRANSA=1 
807:          M=NINTC 
808:          N=1 
809:          K=NCART 
810:          ALPHA=1.0D0 
811:          BETA=0.0D0 
812:          DO I1=1,9 
813:             DESCRA(I1)=0 
814:          ENDDO 
815:          DESCRA(4)=1 
816:          LWORK=NINTC 
817:          IF (PTEST) THEN 
818:             PRINT*,'before dcoomm 1 NEWX:' 
819:             WRITE(*,'(6G20.10)') (NEWX(J1),J1=1,NINTC) 
820:             PRINT*,'before dcoomm 1 CARTX:' 
821:             WRITE(*,'(6G20.10)') (CARTX(J1),J1=1,NCART) 
822:          ENDIF 
823:          CALL DCOOMM(TRANSA,M,N,K,ALPHA,DESCRA,VAL,INDX,JNDX,NNZ, & 
824:      &        NEWX,NINTC,BETA,CARTX,NCART,WORK,LWORK) 
825:          IF (PTEST) THEN 
826:             PRINT*,'after dcoomm 1 NEWX:' 
827:             WRITE(*,'(6G20.10)') (NEWX(J1),J1=1,NINTC) 
828:             PRINT*,'after dcoomm 1 CARTX:' 
829:             WRITE(*,'(6G20.10)') (CARTX(J1),J1=1,NCART) 
830:          ENDIF          
831: !     
832: !    Cholesky decompose GCS 
833: !    use lapack routine DPBTRF 
834: !    LDGCS is leading dimension of GCS 
835: !     
836:          UPLO='U' 
837:          N=NCART 
838:          LDGCS=KD+1 
839:         IF (PTEST) THEN 
840:            PRINT*,'GCS before DPBTRF:' 
841:            WRITE(*,'(6G18.10)') ((GCS(J1,J2),J1=1,KD+1),J2=1,NCART) 
842:         ENDIF 
843: !    CALL MYCPU_TIME(TESTTIME1, .FALSE.) 
844:          CALL DPBTRF(UPLO,N,KD,GCS,LDGCS,INFO) 
845: !    CALL MYCPU_TIME(TESTTIME2, .FALSE.) 
846: !    TOTCHOLTIME = TOTCHOLTIME + TESTTIME2 - TESTTIME1 
847:          IF (INFO.NE.0) PRINT*,'after DPBTRF INFO=',INFO 
848:  
849: !msb50 
850: !        IF (PTEST) THEN 
851: !           PRINT*,'GCS after DPBTRF:' 
852: !           WRITE(*,'(6G18.10)') ((GCS(J1,J2),J1=1,KD+1),J2=1,NCART) 
853: !        ENDIF 
854: !     
855: !    GCS is now the upper triangular matrix after decomposition 
856: !    SOLVE G*DELTAX=CARTX FOR DELTAX, WHERE G IS ORIGINAL BEETBEE+INTEPSILON*I 
857: !    done using lapack routine DPBTRS with input of the upper matrix from 
858: !    cholesky decomposition 
859: !     
860:          NRHS=1 
861:          LDB=NCART 
862:          CALL DPBTRS(UPLO,N,KD,NRHS,GCS,LDGCS,CARTX,LDB,INFO) 
863:          IF (INFO.NE.0) PRINT*,'after DPBTRS INFO=',INFO 
864:         IF (PTEST) THEN 
865:            PRINT*,'GCS after DPBTRS:' 
866:            WRITE(*,'(6G18.10)') ((GCS(J1,J2),J1=1,KD+1),J2=1,NCART) 
867:            PRINT*,'after DPBTRS CARTX:' 
868:            WRITE(*,'(6G20.10)') (CARTX(J1),J1=1,NCART) 
869:         ENDIF 
870:  
871: !     
872: !    CARTX is now solution of CARTX=A*X 
873: !    evaluate XNEW2=B*CARTX 
874: !     
875:  
876:          IF (INTNEWT) THEN 
877:             ! evaluate XNEW2 directly from new cartesian coords 
878:             NOCOOR = .FALSE.; NODERV = .TRUE. 
879:             CALL AMB_GETBEENAT(DUMCOORDS+CARTX,DUMXINT2,DUMVAL,DUMINDX, & 
880:      &           DUMJNDX,NINTC,NCART,NNZ,DUMGCS2,NOCOOR,NODERV,KD) 
881:  
882:             ! align dihedral angles if using primitive internals 
883:             IF (.NOT.NATINT) THEN  
884:                !msb 50 - probably bug for amber 
885:                PRINT*, "Warning - AMB_TRANSBACKDELTA has natint .FALSE.& 
886:                &- are you sure? -> uncomment stop in ambnatinern_extras" 
887:                STOP  
888:                DO I1 = NBOND+NTHETS+1, NINTC 
889:                   IF (DUMXINT2(I1) - QTARGET(I1).GT.PI) THEN 
890:                      DUMXINT2(I1) = DUMXINT2(I1) - TWOPI 
891:                   ELSE IF (DUMXINT2(I1) - QTARGET(I1).LT.-PI) THEN 
892:                      DUMXINT2(I1) = DUMXINT2(I1) + TWOPI 
893:                   ENDIF                   
894:                ENDDO 
895:             ENDIF 
896:             NEWX = QTARGET - DUMXINT2 
897:  
898:          ELSE ! solve linear equation only; iterate away effect of epsilon 
899:             TRANSA=0 
900:             N=1 
901:             CALL DCOOMM(TRANSA,M,N,K,ALPHA,DESCRA,VAL,INDX,JNDX,NNZ, & 
902:      &           CARTX,NCART,BETA,XNEW2,NINTC,WORK,LWORK) 
903:  
904:             IF (PTEST) THEN 
905:                PRINT*,'XNEW2 after DCOOMM 2:' 
906:                WRITE(*,'(6G18.10)') (XNEW2(J1),J1=1,NINTC) 
907:                PRINT*,'after DCOOMM 2 CARTX:' 
908:                WRITE(*,'(6G20.10)') (CARTX(J1),J1=1,NCART) 
909:             ENDIF 
910:             DO I1=1,NINTC 
911:                DUMX(I1)=DUMX(I1)+XNEW2(I1) ! sum of f(x_i) instead of f(sum(x_i)) 
912:                NEWX(I1)=X(I1)-DUMX(I1) 
913:             ENDDO 
914:          ENDIF 
915:  
916: !    project translation out of CARTX 
917:          CARTDOTX=0.0D0 
918:          CARTDOTY=0.0D0 
919:          CARTDOTZ=0.0D0 
920:          DO I1=1,NATOMS 
921:             CARTDOTX=CARTDOTX+CARTX(3*(I1-1)+1) 
922:             CARTDOTY=CARTDOTY+CARTX(3*(I1-1)+2) 
923:             CARTDOTZ=CARTDOTZ+CARTX(3*(I1-1)+3) 
924:          ENDDO 
925:          DO I1=1,NATOMS 
926:             CARTX(3*(I1-1)+1)=CARTX(3*(I1-1)+1)-CARTDOTX/NATOMS 
927:             CARTX(3*(I1-1)+2)=CARTX(3*(I1-1)+2)-CARTDOTY/NATOMS 
928:             CARTX(3*(I1-1)+3)=CARTX(3*(I1-1)+3)-CARTDOTZ/NATOMS 
929:          ENDDO 
930: !    update DUMCOORDS for next call to GETBEE; check for convergence 
931:          SUMDC2=0.0D0 
932:          DO I1=1,NCART 
933:             DUMCOORDS(I1)=DUMCOORDS(I1)+CARTX(I1) 
934:             SUMDC2=SUMDC2+CARTX(I1)*CARTX(I1) 
935:          ENDDO 
936:          MEANDC2=SUMDC2/NINTC 
937:  
938:          IF (MEANDC2.LT.BACKTCUTOFF) THEN                    
939:             LCONVERGE=.TRUE. 
940:             FAILED=.FALSE. 
941:          ENDIF 
942:          NITS=NITS+1 
943:  
944:          IF (PTEST2) PRINT '(A,I5,2G20.10)',& 
945:              'transbackdelta>> NITS,MEANDC2, int erro =',NITS, MEANDC2,& 
946:              SQRT(DOT_PRODUCT(NEWX(1:NINTC),NEWX(1:NINTC))) 
947:  
948:          IF (SQRT(DOT_PRODUCT(NEWX(1:NINTC),NEWX(1:NINTC))).GT.100) THEN 
949:             print*, 'ERROR! deltatransform diverged.' 
950:             CALL DUMPCOORDS(COORDS,'diverged.xyz', .FALSE.) 
951:             OPEN(unit=88,FILE='badstp', STATUS = 'UNKNOWN') 
952:             WRITE(88,'(G25.15)') (X(J1), J1 = 1,NINTC) 
953:             CLOSE(88) 
954:             LCONVERGE = .TRUE. 
955:             FAILED = .TRUE. 
956:             CARTX(:) = 0.0D0 
957:             RETURN 
958:          ENDIF 
959:  
960:          IF (NITS.GT.100) THEN 
961:             PRINT*,'WARNING internal coordinate transformation& 
962:                & did not converge, MEANDC2=',MEANDC2 
963:             FAILED = .TRUE. 
964:             LCONVERGE=.TRUE. 
965:          ENDIF 
966: !     
967:       ENDDO 
968:  
969: !    update CARTX so that it is the difference between final DUMCOORDS and COORDS 
970: !    print *,'tdb here1' 
971:  
972:       DO I1=1,3*NATOMS 
973:          CARTX(I1)=DUMCOORDS(I1)-COORDS(I1) 
974:       ENDDO 
975: !     
976: !    project translation out of CARTX 
977:       CARTDOTX=0.0D0 
978:       CARTDOTY=0.0D0 
979:       CARTDOTZ=0.0D0 
980:       DO I1=1,NATOMS 
981:          CARTDOTX=CARTDOTX+CARTX(3*(I1-1)+1) 
982:          CARTDOTY=CARTDOTY+CARTX(3*(I1-1)+2) 
983:          CARTDOTZ=CARTDOTZ+CARTX(3*(I1-1)+3) 
984:       ENDDO 
985:       DO I1=1,NATOMS 
986:          CARTX(3*(I1-1)+1)=CARTX(3*(I1-1)+1)-CARTDOTX/NATOMS 
987:          CARTX(3*(I1-1)+2)=CARTX(3*(I1-1)+2)-CARTDOTY/NATOMS 
988:          CARTX(3*(I1-1)+3)=CARTX(3*(I1-1)+3)-CARTDOTZ/NATOMS 
989:       ENDDO 
990:  
991:       IF(BACKTPRINTERR) print*, 'transbackdelta>>> error = ', SQRT(DOT_PRODUCT(NEWX(1:NINTC),NEWX(1:NINTC))) 
992:  
993:       RETURN 
994:       END 


r15547/charmmBildc.f 2017-01-21 10:33:12.122250000 +0000 r15546/charmmBildc.f 2017-01-21 10:33:18.966250000 +0000
1149:                       TW_DIFFP(i) = DABS(DIFF_PHI(i))/MAX_DIFF1149:                       TW_DIFFP(i) = DABS(DIFF_PHI(i))/MAX_DIFF
1150:                       NTW_ANGLES = NTW_ANGLES +11150:                       NTW_ANGLES = NTW_ANGLES +1
1151:                    ENDIF1151:                    ENDIF
1152:                    A_OMEGA = A_OMEGA +11152:                    A_OMEGA = A_OMEGA +1
1153:                ENDIF1153:                ENDIF
1154:             ENDIF1154:             ENDIF
1155:          ENDDO1155:          ENDDO
1156:       ENDIF1156:       ENDIF
1157:       PRINT*, "NTW_ANGLES", NTW_ANGLES1157:       PRINT*, "NTW_ANGLES", NTW_ANGLES
1158:       END SUBROUTINE GET_TWISTABLE1158:       END SUBROUTINE GET_TWISTABLE
1159:  
1160:  
1161: !********************************************************************************** 
1162:  
1163:       SUBROUTINE AMB_PATOM(ATOM_NO, ATRES, ATOM) 
1164: !msb50: give it a residue number and e.g. "CD1" and in gives you the number of CD1 
1165:       use modamber9 
1166:       use commons, only : NATOMS 
1167:       IMPLICIT NONE 
1168:       CHARACTER(LEN=4),INTENT(IN) :: ATOM  
1169:       INTEGER,INTENT(IN) :: ATRES 
1170:       INTEGER, INTENT(OUT) :: ATOM_NO 
1171:       INTEGER ::k 
1172:        
1173: !     FOR GOD'S SAKE REMEMBER 
1174: !     ix(i02) STARTS WITH 1 
1175: ! 
1176:       ATOM_NO = -9999 
1177: !     nat_res =ix(i02+ATRES)- ix(i02+ATRES-1)  !number of atoms in current residue 
1178:       DO k=ix(i02+ATRES-1),ix(i02+ATRES) 
1179:           IF (ATOM==ih(m04+k-1)) THEN 
1180:              ATOM_NO = k 
1181:              EXIT 
1182:           ENDIF 
1183:       ENDDO  
1184:       IF (ATOM_NO.LT.0) THEN 
1185:          PRINT*, "ERROR: atom not found", ATOM 
1186:       ENDIF 
1187:       END SUBROUTINE AMB_PATOM 


r15547/charmm_main.src 2017-01-21 10:33:07.575188235 +0000 r15546/charmm_main.src 2017-01-21 10:33:14.886250000 +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
703:         CALL INTCOR(COMLYN,COMLEN)705:         CALL INTCOR(COMLYN,COMLEN)
704:       ELSE IF (WRD.EQ.'INQU') THEN706:       ELSE IF (WRD.EQ.'INQU') THEN
705: C---- Procedure PROCESS-INQUIRE-FILE-COMMAND707: C---- Procedure PROCESS-INQUIRE-FILE-COMMAND
706: C708: C
707: C     inquire all open files709: C     inquire all open files
708:         IF(PRNLEV.GE.2) WRITE(OUTU,'(A)')710:         IF(PRNLEV.GE.2) WRITE(OUTU,'(A)')
709:      $        ' CHARMM: list of open files:'711:      $        ' CHARMM: list of open files:'
710:         DO I=1,99712:         DO I=1,99
711:           CALL VINQRE('UNIT',TNAME,TMAX,TLEN,QOPEN,QFORM,QWRITE,I)713:           CALL VINQRE('UNIT',TNAME,TMAX,TLEN,QOPEN,QFORM,QWRITE,I)
712:           IF (QOPEN) THEN714:           IF (QOPEN) THEN
829:       ELSE IF (WRD.EQ.'POLA') THEN831:       ELSE IF (WRD.EQ.'POLA') THEN
830:         CALL POLAR0832:         CALL POLAR0
831: ##ENDIF833: ##ENDIF
832: ##IF PATHINT834: ##IF PATHINT
833:       ELSE IF (WRD.EQ.'PINT') THEN835:       ELSE IF (WRD.EQ.'PINT') THEN
834:         CALL PINT(COMLYN,COMLEN)836:         CALL PINT(COMLYN,COMLEN)
835: ##ENDIF837: ##ENDIF
836:       ELSE IF (WRD.EQ.'PRES') THEN838:       ELSE IF (WRD.EQ.'PRES') THEN
837:         CALL GETPRS(COMLYN,COMLEN)839:         CALL GETPRS(COMLYN,COMLEN)
838:       ELSE IF (WRD.EQ.'PRIN') THEN840:       ELSE IF (WRD.EQ.'PRIN') THEN
 841:         PRINT*, "in chsetup msb50 before mainio l827"
839:         CALL MAINIO(WRD)842:         CALL MAINIO(WRD)
840:       ELSE IF (WRD.EQ.'PULL') THEN843:       ELSE IF (WRD.EQ.'PULL') THEN
841: C---- 23-Jul-96, LN844: C---- 23-Jul-96, LN
842:         CALL PULL(COMLYN,COMLEN)845:         CALL PULL(COMLYN,COMLEN)
843:       ELSE IF (WRD.EQ.'QUAN') THEN846:       ELSE IF (WRD.EQ.'QUAN') THEN
844:         CALL QMDEFN(COMLYN,COMLEN)847:         CALL QMDEFN(COMLYN,COMLEN)
845:       ELSE IF (WRD.EQ.'READ') THEN848:       ELSE IF (WRD.EQ.'READ') THEN
846:         CALL MAINIO(WRD)849:         CALL MAINIO(WRD)
847: ##IF GENETIC850: ##IF GENETIC
848: C---- 12-Jan-98, CLBIII851: C---- 12-Jan-98, CLBIII


r15547/intcoords.f90 2017-01-21 10:33:12.442250000 +0000 r15546/intcoords.f90 2017-01-21 10:33:19.250250000 +0000
  8:   IMPLICIT NONE  8:   IMPLICIT NONE
  9:   9: 
 10: CONTAINS 10: CONTAINS
 11:     SUBROUTINE INTINTERPOLATE(RS,RF,NIMI,NIMC, XINTERP, PTEST,OUTFAILED) 11:     SUBROUTINE INTINTERPOLATE(RS,RF,NIMI,NIMC, XINTERP, PTEST,OUTFAILED)
 12:     ! do a linear interpolation between start and end  12:     ! do a linear interpolation between start and end 
 13:     ! using internal coordinates, with NIMI images 13:     ! using internal coordinates, with NIMI images
 14:     ! then place NIMC images at equidistant cartesian 14:     ! then place NIMC images at equidistant cartesian
 15:     ! points along the resulting interpolated path     15:     ! points along the resulting interpolated path    
 16:     ! results go in XCART 16:     ! results go in XCART
 17:  17: 
 18:     USE KEY, ONLY : RIGIDBODY, TWOD, BULKT, AMBERT 18:     USE KEY, ONLY : RIGIDBODY, TWOD, BULKT
 19:  19: 
 20:     IMPLICIT NONE 20:     IMPLICIT NONE
 21:  21: 
 22:     DOUBLE PRECISION, INTENT(IN) :: RS(3*NATOMS), RF(3*NATOMS) 22:     DOUBLE PRECISION, INTENT(IN) :: RS(3*NATOMS), RF(3*NATOMS)
 23:     INTEGER, INTENT(IN) :: NIMC 23:     INTEGER, INTENT(IN) :: NIMC
 24:     INTEGER, INTENT(INOUT) :: NIMI 24:     INTEGER, INTENT(INOUT) :: NIMI
 25:     DOUBLE PRECISION, INTENT(OUT) :: XINTERP(3*NATOMS*NIMC) 25:     DOUBLE PRECISION, INTENT(OUT) :: XINTERP(3*NATOMS*NIMC)
 26:     LOGICAL, INTENT(IN) :: PTEST 26:     LOGICAL, INTENT(IN) :: PTEST
 27:     LOGICAL, INTENT(OUT) :: OUTFAILED 27:     LOGICAL, INTENT(OUT) :: OUTFAILED
 28:  28: 
 35:     DOUBLE PRECISION :: EMAX, EMAXC, ENERGY, GDUMMY(3*NATOMS), TMPRMS 35:     DOUBLE PRECISION :: EMAX, EMAXC, ENERGY, GDUMMY(3*NATOMS), TMPRMS
 36:     DOUBLE PRECISION :: XCART2(3*NATOMS*NIMC), CARTDIFF(3*NATOMS), GINTDUMMY(NINTC) 36:     DOUBLE PRECISION :: XCART2(3*NATOMS*NIMC), CARTDIFF(3*NATOMS), GINTDUMMY(NINTC)
 37:     INTEGER :: ISTAT 37:     INTEGER :: ISTAT
 38:     LOGICAL :: FAILED 38:     LOGICAL :: FAILED
 39:     DOUBLE PRECISION :: RFORIG(3*NATOMS) 39:     DOUBLE PRECISION :: RFORIG(3*NATOMS)
 40:  40: 
 41:     ! alignment stuff 41:     ! alignment stuff
 42:     DOUBLE PRECISION :: DISTF, DIST, DIST2, RMAT(3,3) 42:     DOUBLE PRECISION :: DISTF, DIST, DIST2, RMAT(3,3)
 43:     CHARACTER(LEN=5) :: ZSYMSAVE 43:     CHARACTER(LEN=5) :: ZSYMSAVE
 44:     COMMON /SYS/ ZSYMSAVE 44:     COMMON /SYS/ ZSYMSAVE
 45:     DOUBLE PRECISION TEST(3*NATOMS) 
 46:  45: 
 47:     IF (.NOT.NATINT) THEN 46:     IF (.NOT.NATINT) THEN
 48:        print*, 'intinterpolate >> INTINTERPOLATE is not set up to work without NATINT' 47:        print*, 'intinterpolate >> INTINTERPOLATE is not set up to work without NATINT'
 49:        STOP 48:        STOP
 50:     ENDIF     49:     ENDIF    
 51:  50: 
 52:     OUTFAILED = .FALSE. 51:     OUTFAILED = .FALSE.
 53:  52: 
 54:     BACKTCUTOFF = INTERPBACKTCUT 53:     BACKTCUTOFF = INTERPBACKTCUT
 55:  54: 
 72:           CALL POTENTIAL(XCART2(NC*(IM-1)+1:NC*IM), ENERGY, GDUMMY, .FALSE., .FALSE., & 71:           CALL POTENTIAL(XCART2(NC*(IM-1)+1:NC*IM), ENERGY, GDUMMY, .FALSE., .FALSE., &
 73:                & TMPRMS, .FALSE., .FALSE.) 72:                & TMPRMS, .FALSE., .FALSE.)
 74:            73:           
 75:           IF (ENERGY.GT.EMAXC) EMAXC = ENERGY 74:           IF (ENERGY.GT.EMAXC) EMAXC = ENERGY
 76:        ENDDO 75:        ENDDO
 77:     ENDIF 76:     ENDIF
 78:  77: 
 79:     IF (NATINT) PREVDIH => DIHINFO(1,:) 78:     IF (NATINT) PREVDIH => DIHINFO(1,:)
 80:     ALIGNDIR = .FALSE. 79:     ALIGNDIR = .FALSE.
 81:     CALL CART2INT(RS,RSINT) 80:     CALL CART2INT(RS,RSINT)
 82:     PRINT*, "intcoords> start coor" 81: 
 83: !   PRINT*, RSINT(106:310) 
 84:     
 85:     ! align individual dihedrals of all other images and endpoint to start 82:     ! align individual dihedrals of all other images and endpoint to start
 86:     DO IM = 2,NIMI+2 83:     DO IM = 2,NIMI+2
 87:        DIHINFO(IM,:) = DIHINFO(1,:) 84:        DIHINFO(IM,:) = DIHINFO(1,:)
 88:     ENDDO 85:     ENDDO
 89:     PREVDIH => DIHINFO(NIMI+2,:) 86:     PREVDIH => DIHINFO(NIMI+2,:)
 90:     ALIGNDIR = .TRUE. 87:     ALIGNDIR = .TRUE.
 91:     CALL CART2INT(RF,RFINT) 88:     CALL CART2INT(RF,RFINT)
 92:     PRINT*, "intcoords> int coords  fin coord" 89:     
 93: !    PRINT*, RFINT(106) 
 94:  
 95:     ALIGNDIR = .FALSE. 90:     ALIGNDIR = .FALSE.
 96:      91:     
 97:     INTDIFF = RFINT(:) - RSINT(:) 92:     INTDIFF = RFINT(:) - RSINT(:)
 98:  93: 
 99:     ! place the images in internals and convert to cartesians 94:     ! place the images in internals and convert to cartesians
100:     IF (PTEST) THEN 95:     IF (PTEST) THEN
101:        print*, 'intinterpolate>> step size: ', SQRT(DOT_PRODUCT(INTDIFF,INTDIFF)) / (NIMI+1) 96:        print*, 'intinterpolate>> step size: ', SQRT(DOT_PRODUCT(INTDIFF,INTDIFF)) / (NIMI+1)
102: !      CALL FLUSH(6,ISTAT) 97: !      CALL FLUSH(6,ISTAT)
103:     ENDIF 98:     ENDIF
104:  99: 
105:     XCART(0,1:NC) = RS(:)100:     XCART(0,1:NC) = RS(:)
106:     ARCS(0) = 0.0D0101:     ARCS(0) = 0.0D0
107:     DO IM = 1,NIMI102:     DO IM = 1,NIMI
108:        PREVDIH => DIHINFO(IM+1,:)103:        PREVDIH => DIHINFO(IM+1,:)
109:        IF (AMBERT) THEN !msb50 
110:           CALL AMB_TRANSBACKDELTA(INTDIFF(:) * 1.0D0 / (NIMI+1),DIFFS(IM,:), & 
111:    &         XCART(IM-1,:),NINTC,3*NATOMS,NNZ,KD,FAILED,.FALSE.,INTEPSILON) 
112:        ELSE 
113:        CALL TRANSBACKDELTA(INTDIFF(:) * 1.0D0 / (NIMI+1),DIFFS(IM,:),XCART(IM-1,:),&104:        CALL TRANSBACKDELTA(INTDIFF(:) * 1.0D0 / (NIMI+1),DIFFS(IM,:),XCART(IM-1,:),&
114:             & NINTC,3*NATOMS,NNZ,KD,FAILED,.FALSE.,INTEPSILON)       105:             & NINTC,3*NATOMS,NNZ,KD,FAILED,.FALSE.,INTEPSILON)          
115:        ENDIF 
116:        XCART(IM,:) = XCART(IM-1,:) + DIFFS(IM,:)                        106:        XCART(IM,:) = XCART(IM-1,:) + DIFFS(IM,:)                        
117: 107: 
118:        ! line up structure with start point108:        ! line up structure with start point
119:        CALL NEWMINDIST(RS(:),XCART(IM,:),NATOMS,DIST,.FALSE.,.FALSE.,ZSYM(1),.FALSE.,.FALSE.,.FALSE.,RMAT)109:        CALL NEWMINDIST(RS(:),XCART(IM,:),NATOMS,DIST,.FALSE.,.FALSE.,ZSYM(1),.FALSE.,.FALSE.,.FALSE.,RMAT)
120: 110: 
121:        DIFFS(IM,:) = XCART(IM,:) - XCART(IM-1,:)111:        DIFFS(IM,:) = XCART(IM,:) - XCART(IM-1,:)
122:        DISTS(IM) = SQRT(DOT_PRODUCT(DIFFS(IM,:), DIFFS(IM,:)))112:        DISTS(IM) = SQRT(DOT_PRODUCT(DIFFS(IM,:), DIFFS(IM,:)))
123:        ARCS(IM) = ARCS(IM-1) + DISTS(IM)113:        ARCS(IM) = ARCS(IM-1) + DISTS(IM)
124: 114: 
125:        IF (INTERPSIMPLE) XINTERP(NC*(IM-1)+1:NC*IM) = XCART(IM,:)115:        IF (INTERPSIMPLE) XINTERP(NC*(IM-1)+1:NC*IM) = XCART(IM,:)
190:     IF (INTERPCHOICE) THEN180:     IF (INTERPCHOICE) THEN
191:        IF (EMAX.GT.EMAXC) THEN181:        IF (EMAX.GT.EMAXC) THEN
192:           IF (PTEST) print '(A,2G20.10)', 'intinterpolate>> using cartesian interpolation', EMAXC, EMAX182:           IF (PTEST) print '(A,2G20.10)', 'intinterpolate>> using cartesian interpolation', EMAXC, EMAX
193:           XINTERP(:) = XCART2(:)183:           XINTERP(:) = XCART2(:)
194:        ELSE184:        ELSE
195:           IF (PTEST) print '(A,2G20.10)', 'intinterpolate>> using internal interpolation', EMAXC, EMAX185:           IF (PTEST) print '(A,2G20.10)', 'intinterpolate>> using internal interpolation', EMAXC, EMAX
196:        ENDIF186:        ENDIF
197:     ENDIF187:     ENDIF
198: !   CALL FLUSH(6,ISTAT)188: !   CALL FLUSH(6,ISTAT)
199: 189: 
200:     PRINT*, "intinterp> deall" 
201:     DEALLOCATE(XCART, DIFFS, DISTS, ARCS, DIHINFO)190:     DEALLOCATE(XCART, DIFFS, DISTS, ARCS, DIHINFO)
202: 191: 
203:     IF (PTEST) print*, 'Finished initial interpolation'192:     IF (PTEST) print*, 'Finished initial interpolation'
204:     RETURN193:     RETURN
205: 194: 
206:   END SUBROUTINE INTINTERPOLATE195:   END SUBROUTINE INTINTERPOLATE
207: 196: 
208:   SUBROUTINE INTMINPERM(RS,RF,DISTANCE,RMAT,PTEST)197:   SUBROUTINE INTMINPERM(RS,RF,DISTANCE,RMAT,PTEST)
209:     ! minimise Cartesian distance without permuting then198:     ! minimise Cartesian distance without permuting then
210:     ! permute the permutable groups to minimize distance in single torsions199:     ! permute the permutable groups to minimize distance in single torsions
211:     ! changes RF to resulting best alignment with RS200:     ! changes RF to resulting best alignment with RS
212:     ! output distance is cartesian distance **squared**201:     ! output distance is cartesian distance **squared**
213: 202: 
214:     USE KEY, ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSWAP, SWAP1, SWAP2, AMBERT203:     USE KEY, ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSWAP, SWAP1, SWAP2
215:     USE MODAMBER9, ONLY : ih, ix, NRES, i02,m04,m02 
216: 204: 
217:     IMPLICIT NONE205:     IMPLICIT NONE
218: 206: 
219:     DOUBLE PRECISION, INTENT(IN) :: RS(3*NATOMS)207:     DOUBLE PRECISION, INTENT(IN) :: RS(3*NATOMS)
220:     DOUBLE PRECISION, INTENT(INOUT) :: RF(3*NATOMS), RMAT(3,3)208:     DOUBLE PRECISION, INTENT(INOUT) :: RF(3*NATOMS), RMAT(3,3)
221:     DOUBLE PRECISION, INTENT(OUT) :: DISTANCE209:     DOUBLE PRECISION, INTENT(OUT) :: DISTANCE
222:     LOGICAL, INTENT(IN) :: PTEST210:     LOGICAL, INTENT(IN) :: PTEST
223: 211: 
224:     INTEGER :: PG, SW, I, J, K, A1, A2, A3, B1, B2, B3, START212:     INTEGER :: PG, SW, I, J, K, A1, A2, A3, B1, B2, B3, START
225:     DOUBLE PRECISION :: RSINT(NINTC), RFINT(NINTC), INTDIFF(NDIH)213:     DOUBLE PRECISION :: RSINT(NINTC), RFINT(NINTC), INTDIFF(NDIH)
226:     DOUBLE PRECISION :: INTDIST, CURINTDIST214:     DOUBLE PRECISION :: INTDIST, CURINTDIST
227:     DOUBLE PRECISION :: TMP(3), TMP2(3), TMP3(3), STARTDIH(NDIH), RFTMP(3*NATOMS)215:     DOUBLE PRECISION :: TMP(3), TMP2(3), TMP3(3), STARTDIH(NDIH), RFTMP(3*NATOMS)
228:     DOUBLE PRECISION :: RFTMP2(3*NATOMS-6)216: 
229:     ! alignment stuff217:     ! alignment stuff
 218:     DOUBLE PRECISION :: DISTF
230:     CHARACTER(LEN=5) :: ZSYMSAVE219:     CHARACTER(LEN=5) :: ZSYMSAVE
231:     COMMON /SYS/ ZSYMSAVE220:     COMMON /SYS/ ZSYMSAVE
232:     DOUBLE PRECISION :: DISTANCE2 
233:     INTEGER :: III 
234:     INTEGER ::glyc 
235:     LOGICAL ::glyswap 
236:  
237:     glyc = 0 
238:     glyswap = .FALSE. 
239:     IF (AMBERT) THEN !msb50 - for amber9.ff03 and glycine this doesn't work 
240:       DO III = 1,NRES!as swapping the hydrogens makes no difference in natural internals in this case 
241:          IF (ih(m02+III-1).EQ.'GLY'.OR.ih(m02+III-1).EQ.'NGLY'.OR.ih(m02+III-1).EQ.'CGLY') THEN 
242:              glyc = III 
243:              PRINT*, "intminperm> I hope you know what you are doing" 
244:          ENDIF 
245:          EXIT 
246:       ENDDO 
247: !      PRINT*, "intminperm> May not be the smartest idea with Amberff03 unless you are sure it is" 
248:     ENDIF 
249:  
250:     DISTANCE2 = DOT_PRODUCT(RF-RS,RF-RS) 
251:     PRINT*, "intminperm> initial cart dist", DISTANCE2 
252: 221: 
253:     ! align without permuting222:     ! align without permuting
254:     CALL NEWMINDIST(RS(:),RF(:),NATOMS,DISTANCE,.FALSE.,.FALSE.,ZSYM(1),.FALSE.,.FALSE.,.FALSE.,RMAT)223:     CALL NEWMINDIST(RS(:),RF(:),NATOMS,DISTANCE,.FALSE.,.FALSE.,ZSYM(1),.FALSE.,.FALSE.,.FALSE.,RMAT)
255: 224: 
256:     PREVDIH => DIHINFOSINGLE(:)225:     PREVDIH => DIHINFOSINGLE(:)
257: 226: 
258:     PREVDIH(:) = 0.0D0227:     PREVDIH(:) = 0.0D0
259:     CALL CART2INT(RS, RSINT)228:     CALL CART2INT(RS, RSINT)
260: 229: 
261:     STARTDIH(:) = PREVDIH(:)230:     STARTDIH(:) = PREVDIH(:)
262:     ALIGNDIR = .FALSE.231:     ALIGNDIR = .FALSE.
263: 232: 
264:     CALL CART2INT(RF, RFINT)    233:     CALL CART2INT(RF, RFINT)    
265:     RFTMP2(:) = RFINT(:)  
266: 234: 
267:     INTDIFF(:) = PREVDIH(:)-STARTDIH(:)235:     INTDIFF(:) = PREVDIH(:)-STARTDIH(:)
268:     CURINTDIST = SQRT(DOT_PRODUCT(INTDIFF, INTDIFF))236:     CURINTDIST = SQRT(DOT_PRODUCT(INTDIFF, INTDIFF))
269:     IF (PTEST) print*, 'Starting intdistance: ', CURINTDIST237:     IF (PTEST) print*, 'Starting intdistance: ', CURINTDIST
270: 238: 
271:     START = 1239:     START = 1
272:     DO PG = 1,NPERMGROUP240:     DO PG = 1,NPERMGROUP
273:  241: 
274:        IF (NPERMSIZE(PG).EQ.2) THEN242:        IF (NPERMSIZE(PG).EQ.2) THEN
275:           RFTMP(:) = RF(:)243:           RFTMP(:) = RF(:)
276: 244: 
277:           ! A1 and A2 are the atom numbers to permute245:           ! A1 and A2 are the atom numbers to permute
278:           A1 = PERMGROUP(START)246:           A1 = PERMGROUP(START)
279:           A2 = PERMGROUP(START+1)247:           A2 = PERMGROUP(START+1)
280: 248: 
281:           IF (PTEST) PRINT*, "which PERMGROUP", PG 
282:           IF (PTEST) print*, 'Swapping atoms: ', A1, A2249:           IF (PTEST) print*, 'Swapping atoms: ', A1, A2
283: 250: 
284:           IF (glyc.NE.0) THEN 
285:               glyswap = .FALSE. 
286:               IF (ix(i02+glyc-1).LE.A1 .AND.(ix(i02+glyc).GT.A1).AND.& 
287:      &           (ih(m04+A1-1).EQ.'HA2 '.OR.ih(m04+A1-1).EQ.'HA3 ')) THEN 
288:                 !the second atom has to be A2 in gly 
289:                 IF (ih(m04+A2-1).EQ.'HA3 '.OR. ih(m04+A2-1).EQ.'HA2 ') THEN 
290:                    glyswap = .TRUE. 
291:                    RF(3*(A1-1)+1:3*A1) = RFTMP(3*(A2-1)+1:3*A2) 
292:                    RF(3*(A2-1)+1:3*A2) = RFTMP(3*(A1-1)+1:3*A1) 
293:                    !move any other groups that must be dragged along 
294:                    IF (NSWAP(PG).NE.0) THEN 
295:                       PRINT*, "intminperm>this is not glycine. you DON'T know&  
296:      &                     what you're doing" 
297:                       STOP 
298:                    ENDIF 
299:                    DISTANCE = DOT_PRODUCT(RFTMP-RS,RFTMP-RS) 
300:                    DISTANCE2 = DOT_PRODUCT(RF-RS,RF-RS) 
301:                    IF (DISTANCE2.LT.DISTANCE) THEN 
302:                     IF (PTEST) print '(a40, f17.5,a10,f17.5)',"keep permutation& 
303:      & according to cartesians (gly): dist1",DISTANCE, "dist_new", DISTANCE2 
304:                    ELSE 
305:                       RF(:) = RFTMP(:) !unddo 
306:                    ENDIF 
307:                  ELSE 
308:                     PRINT*, "intminperm> sth wrong with GLY" 
309:                     STOP 
310:                  ENDIF 
311:               ENDIF 
312:               IF (glyswap .EQV. .TRUE.) THEN 
313:                 START = START + NPERMSIZE(PG) 
314:                 CYCLE 
315:                 ENDIF 
316:           ENDIF 
317:  
318:           RF(3*(A1-1)+1:3*A1) = RFTMP(3*(A2-1)+1:3*A2)251:           RF(3*(A1-1)+1:3*A1) = RFTMP(3*(A2-1)+1:3*A2)
319:           RF(3*(A2-1)+1:3*A2) = RFTMP(3*(A1-1)+1:3*A1)252:           RF(3*(A2-1)+1:3*A2) = RFTMP(3*(A1-1)+1:3*A1)
320: 253: 
321:           !move any other groups that must be dragged along254:           !move any other groups that must be dragged along
322:           DO SW = 1,NSWAP(PG)255:           DO SW = 1,NSWAP(PG)
323:              A1 = SWAP1(PG,SW)256:              A1 = SWAP1(PG,SW)
324:              A2 = SWAP2(PG,SW)257:              A2 = SWAP2(PG,SW)
325: 258: 
326:              IF (PTEST) print*, 'Dragging swap: ', A1, A2259:              IF (PTEST) print*, 'Dragging swap: ', A1, A2
327: 260: 
330:           ENDDO263:           ENDDO
331: 264: 
332:           PREVDIH(:) = STARTDIH(:)265:           PREVDIH(:) = STARTDIH(:)
333:           CALL CART2INT(RF,RFINT)266:           CALL CART2INT(RF,RFINT)
334: 267: 
335:           INTDIFF(:) = PREVDIH(:) - STARTDIH(:)268:           INTDIFF(:) = PREVDIH(:) - STARTDIH(:)
336:           INTDIST = SQRT(DOT_PRODUCT(INTDIFF, INTDIFF))269:           INTDIST = SQRT(DOT_PRODUCT(INTDIFF, INTDIFF))
337: 270: 
338:           IF (PTEST) print*, 'torsion distance: ', INTDIST271:           IF (PTEST) print*, 'torsion distance: ', INTDIST
339: 272: 
340:           IF (INTDIST.GE.CURINTDIST) THEN273:           IF (INTDIST.GT.CURINTDIST) THEN
341:              ! undo permutation274:              ! undo permutation
342:              RF(:) = RFTMP(:)275:              RF(:) = RFTMP(:)
343:           ELSE276:           ELSE
344:              IF (PTEST) print*, 'keep permutation'277:              IF (PTEST) print*, 'keep permutation'
345:              CURINTDIST = INTDIST278:              CURINTDIST = INTDIST
346:           ENDIF279:           ENDIF
347:        ELSE IF (NPERMSIZE(PG).EQ.3) THEN280:        ELSE IF (NPERMSIZE(PG).EQ.3) THEN
348:           ! this is an inefficient way of listing permutations281:           ! this is an inefficient way of listing permutations
349:           ! and should be fixed and generalized at some point282:           ! and should be fixed and generalized at some point
350:           A1 = PERMGROUP(START)283:           A1 = PERMGROUP(START)
402: 335: 
403:     DISTANCE = DOT_PRODUCT(RF-RS,RF-RS)336:     DISTANCE = DOT_PRODUCT(RF-RS,RF-RS)
404: 337: 
405:     RETURN338:     RETURN
406:   END SUBROUTINE INTMINPERM339:   END SUBROUTINE INTMINPERM
407: 340: 
408:   SUBROUTINE CART2INT(XCART, XINT) 341:   SUBROUTINE CART2INT(XCART, XINT) 
409:     ! converts cartesians to internals342:     ! converts cartesians to internals
410:     ! using parameters from intcommons343:     ! using parameters from intcommons
411:     ! no derivatives involved344:     ! no derivatives involved
412:     USE KEY, ONLY:AMBERT 
413: 345: 
414:     IMPLICIT NONE346:     IMPLICIT NONE
415:     DOUBLE PRECISION, INTENT(IN) :: XCART(3*NATOMS)347:     DOUBLE PRECISION, INTENT(IN) :: XCART(3*NATOMS)
416:     DOUBLE PRECISION, INTENT(OUT) :: XINT(NINTC)348:     DOUBLE PRECISION, INTENT(OUT) :: XINT(NINTC)
417:     DOUBLE PRECISION :: GDUMMY(3*NATOMS), GINT(NINTC)349:     DOUBLE PRECISION :: GDUMMY(3*NATOMS), GINT(NINTC)
418: 350: 
419:     IF (AMBERT) THEN351:     CALL TRANSFORM(XCART,GDUMMY,XINT,GINT,NINTC,3*NATOMS,NNZ,.FALSE.,.TRUE.,KD,INTEPSILON)
420:         CALL AMBTRANSFORM(XCART,GDUMMY,XINT,GINT,NINTC,3*NATOMS,NNZ,.FALSE.,.TRUE.,KD,INTEPSILON) 
421:     ELSE 
422:         CALL TRANSFORM(XCART,GDUMMY,XINT,GINT,NINTC,3*NATOMS,NNZ,.FALSE.,.TRUE.,KD,INTEPSILON) 
423:     ENDIF 
424: 352: 
425:     RETURN353:     RETURN
426:   END SUBROUTINE CART2INT354:   END SUBROUTINE CART2INT
427: 355: 
428:   SUBROUTINE INTSETUP356:   SUBROUTINE INTSETUP
429:     USE MODAMBER9 
430:     USE KEY, ONLY:AMBERT, DEBUG 
431:     IMPLICIT NONE357:     IMPLICIT NONE
432:     ! set up some of the initial global arrays and constants for working with internals358:     ! set up some of the initial global arrays and constants for working with internals
433:     ! interface to CHARMM    359:     ! interface to CHARMM    
434: 360: 
435:     ALLOCATE(USECART(NATOMS))361:     ALLOCATE(USECART(NATOMS))
436:     USECART(:) = .FALSE.362:     USECART(:) = .FALSE.
437: 363: 
438:     IF (CHRMMT) THEN364:     IF (CHRMMT) THEN
439:        CALL GETCHRMINTPARAM365:        CALL GETCHRMINTPARAM
440:        ALLOCATE(BACKBONE(NATOMS))366:        ALLOCATE(BACKBONE(NATOMS))
441:        CALL GETBACKBONE367:        CALL GETBACKBONE
442:     ENDIF368:     ENDIF
443: 369: 
444:     ! NINTS is the number of internals in commons file; NINTC should be same number in intcommons file370:     ! NINTS is the number of internals in commons file; NINTC should be same number in intcommons file
445:     IF (NATINT) THEN371:     IF (NATINT) THEN
446:        IF (AMBERT) THEN372:        CALL NATINTSETUP
447:           CALL AMB_NATINTSETUP 
448:        ELSE  
449:           CALL NATINTSETUP 
450:        ENDIF 
451:        CALL SETCARTATMS373:        CALL SETCARTATMS
452:        IF (USEPARFILE) THEN374:        IF (USEPARFILE) THEN
453:           CALL GETNATINTERNFILE375:           CALL GETNATINTERNFILE
454:        ELSE376:        ELSE
455:           IF (AMBERT) THEN377:           CALL GETNATINTERN
456:              CALL AMBGETNATINTERN 
457:              IF (DEBUG) PRINT*, "intcoords> after ambgetnatintern" 
458:           ELSE 
459:               CALL GETNATINTERN 
460:           ENDIF 
461:        ENDIF378:        ENDIF
462:        NINTS = NINTC379:        NINTS = NINTC
463: 380: 
464:        ALLOCATE(DIHINFOSINGLE(NDIH))381:        ALLOCATE(DIHINFOSINGLE(NDIH))
465:        ALIGNDIR = .FALSE.382:        ALIGNDIR = .FALSE.
466:        DIHINFOSINGLE(:) = 0.0D0383:        DIHINFOSINGLE(:) = 0.0D0
467:        PREVDIH => DIHINFOSINGLE(:)384:        PREVDIH => DIHINFOSINGLE(:)
468:     ENDIF385:     ENDIF
469: 386: 
470:     IF( AMBERT) THEN387:     CALL GETNNZ(NNZ)
471:        IF (NATINT) THEN 388:     CALL GETKD(KD)
472:           CALL AMBGETNNZNAT(NNZ) 
473:           CALL AMB_GETKDNAT(KD) 
474:        ELSE 
475:           PRINT*, "error in intcoords - natint not set for AMBER" 
476:           STOP 
477:        ENDIF 
478:     ELSE 
479:        CALL GETNNZ(NNZ) 
480:        CALL GETKD(KD) 
481:     ENDIF 
482:     BACKTCUTOFF = MINBACKTCUT389:     BACKTCUTOFF = MINBACKTCUT
483: 390: 
484:     CALL KEYINTPRINT391:     CALL KEYINTPRINT
485:   END SUBROUTINE INTSETUP392:   END SUBROUTINE INTSETUP
486:   393:   
487:   SUBROUTINE SETCARTATMS394:   SUBROUTINE SETCARTATMS
488:     ! pick out the atoms for which to use cartesian coordinates395:     ! pick out the atoms for which to use cartesian coordinates
489:     ! make sure NATINTSETUP is done before this!396:     ! make sure NATINTSETUP is done before this!
490: 397: 
491:     IMPLICIT NONE398:     IMPLICIT NONE
492:     INTEGER :: A399:     INTEGER :: A
493: 400: 
494:     NCRT = 0401:     NCRT = 0
495:     IF (.NOT. ALLOCATED(CARTATMS)) ALLOCATE(CARTATMS(NATOMS))402:     ALLOCATE(CARTATMS(NATOMS))
496: 403: 
497:     DO A = 1,NATOMS404:     DO A = 1,NATOMS
498:        IF (A.GE.CARTATMSTART.OR.(CHRMMT.AND.BBCART.AND.BACKBONE(A))) THEN405:        IF (A.GE.CARTATMSTART.OR.(CHRMMT.AND.BBCART.AND.BACKBONE(A))) THEN
499:           NCRT = NCRT + 1406:           NCRT = NCRT + 1
500:           CARTATMS(NCRT) = A407:           CARTATMS(NCRT) = A
501:           USECART(A) = .TRUE.408:           USECART(A) = .TRUE.
502:        ENDIF       409:        ENDIF       
503:     ENDDO410:     ENDDO
504: 411: 
505:   END SUBROUTINE SETCARTATMS412:   END SUBROUTINE SETCARTATMS
506: 413: 
507:   SUBROUTINE NATINTSETUP414:   SUBROUTINE NATINTSETUP
508:     ! setup up various arrays necessary for natural internals415:     ! setup up various arrays necessary for natural internals
509:     ! this is the routine that defines the natural internal coordinates in terms of angle and torsion coefficients416:     ! this is the routine that defines the natural internal coordinates in terms of angle and torsion coefficients
510: 417: 
511:     IMPLICIT NONE418:     IMPLICIT NONE
512:     INTEGER, PARAMETER :: LARGENEG = -99999999419:     INTEGER, PARAMETER :: LARGENEG = -99999999
513:     DOUBLE PRECISION :: S6R, S2R, S26R, S18R, S12R, A, B, SABR, SABR2420:     DOUBLE PRECISION :: S6R, S2R, S26R, S18R, S12R, A, B, SABR, SABR2
514: 421: 
515:  
516:     PRINT*, "NATOMS", NATOMS 
517:     NBDS = 0; NCNT = 0; NRNG = 0; NFRG = 0; NLDH = 0; NIMP = 0422:     NBDS = 0; NCNT = 0; NRNG = 0; NFRG = 0; NLDH = 0; NIMP = 0
518:     NCRT = 0; NDIH = 0423:     NCRT = 0; NDIH = 0
519: 424: 
520:     ALLOCATE(CENTERS(NATOMS,0:5), RINGS(TOTRES,6), FRINGS(TOTRES,6), LINDIH(NATOMS,-1:8), IMPDIH(NATOMS,4), CARTATMS(NATOMS))425:     ALLOCATE(CENTERS(NATOMS,0:5), RINGS(TOTRES,6), FRINGS(TOTRES,6), LINDIH(NATOMS,-1:8), IMPDIH(NATOMS,4), CARTATMS(NATOMS))
521: 426: 
522:     IF (CARTRESSTART.GT.TOTRES+1.OR.CARTRESSTART.LT.0) THEN427:     IF (CARTRESSTART.GT.TOTRES+1.OR.CARTRESSTART.LT.0) THEN
523:        CARTRESSTART = TOTRES + 1428:        CARTRESSTART = TOTRES + 1
524:        CARTATMSTART = NATOMS + 1429:        CARTATMSTART = NATOMS + 1
525:     ELSE430:     ELSE
526:        CARTATMSTART = RESSTARTS(CARTRESSTART)+1431:        CARTATMSTART = RESSTARTS(CARTRESSTART)+1
554:     COEFF(4,1,1:3)=(/-1,-1,2/)*S6R459:     COEFF(4,1,1:3)=(/-1,-1,2/)*S6R
555:     COEFF(4,2,1:3)=(/1,-1,0/)*S2R460:     COEFF(4,2,1:3)=(/1,-1,0/)*S2R
556: 461: 
557:     COEFF(5,1,1:3)=(/0,1,-1/)*S2R462:     COEFF(5,1,1:3)=(/0,1,-1/)*S2R
558:     COEFF(5,2,1:3)=(/2,-1,-1/)*S6R463:     COEFF(5,2,1:3)=(/2,-1,-1/)*S6R
559: 464: 
560:     COEFF(6,1,1:1) = (/1/)465:     COEFF(6,1,1:1) = (/1/)
561: 466: 
562:     COEFF(7,1,1:2) = (/1,-1/)*S2R467:     COEFF(7,1,1:2) = (/1,-1/)*S2R
563: 468: 
564:     COEFF(8,1,1:5) = (/1,-1,1,-1,0/)*0.5469:     COEFF(8,1,1:3) = (/1,4,1/)*S18R
565:     COEFF(8,2,1:5) = (/-1,-1,1,1,0/)*0.5470:     COEFF(8,2,1:3) = (/1,1,4/)*S18R
566:     COEFF(8,3,1:5) = (/-1,1,1,-1,0/)*0.5471:     COEFF(8,3,1:3) = (/4,1,1/)*S18R
567:     COEFF(8,4,1:5) = (/0,0,0,0,1/) 
568: 472: 
569:     A = COS(144.0*PI/180.0)473:     A = COS(144.0*PI/180.0)
570:     B = COS(72.0*PI/180.0)474:     B = COS(72.0*PI/180.0)
571:     SABR = 1/SQRT(1+2*A*A+2*B*B)475:     SABR = 1/SQRT(1+2*A*A+2*B*B)
572:     SABR2 = 1/SQRT(2*(A-B)**2+2*(1-A)**2)476:     SABR2 = 1/SQRT(2*(A-B)**2+2*(1-A)**2)
573:     COEFF(25,1,1:5) = (/1.0D0,A,B,B,A/)*SABR477:     COEFF(25,1,1:5) = (/1.0D0,A,B,B,A/)*SABR
574:     COEFF(25,2,1:5) = (/0.0D0,A-B,1.0D0-A,A-1.0D0,B-A/)*SABR2478:     COEFF(25,2,1:5) = (/0.0D0,A-B,1.0D0-A,A-1.0D0,B-A/)*SABR2
575:     COEFF(25,3,1:5) = (/B,A,1.0D0,A,B/)*SABR479:     COEFF(25,3,1:5) = (/B,A,1.0D0,A,B/)*SABR
576:     COEFF(25,4,1:5) = (/A-1.0D0,B-A,0.0D0,A-B,1.0D0-A/)*SABR2480:     COEFF(25,4,1:5) = (/A-1.0D0,B-A,0.0D0,A-B,1.0D0-A/)*SABR2
577: 481: 


r15547/intcor2.src 2017-01-21 10:33:08.135188235 +0000 r15546/intcor2.src 2017-01-21 10:33:15.426250000 +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"
151:       IF(.NOT.LFILL) THEN152:       IF(.NOT.LFILL) THEN
152: C153: C
 154:          PRINT*, "hello msb50 in intcor2"
153:          IF(LCLEAR) THEN155:          IF(LCLEAR) THEN
154:             DO I=1,LENIC156:             DO I=1,LENIC
155:                IF(IAR(I).LE.0) IAR(I)=MARK2157:                IF(IAR(I).LE.0) IAR(I)=MARK2
156:                IF(JAR(I).LE.0) IAR(I)=MARK2158:                IF(JAR(I).LE.0) IAR(I)=MARK2
157:                IF(KAR(I).LE.0) IAR(I)=MARK2159:                IF(KAR(I).LE.0) IAR(I)=MARK2
158:                IF(LAR(I).LE.0) IAR(I)=MARK2160:                IF(LAR(I).LE.0) IAR(I)=MARK2
159:             ENDDO161:             ENDDO
160:          ENDIF162:          ENDIF
161: C163: C
162:          II=0164:          II=0
176:                   T1IC(II)=T1IC(I)178:                   T1IC(II)=T1IC(I)
177:                   T2IC(II)=T2IC(I)179:                   T2IC(II)=T2IC(I)
178:                   PIC(II)=PIC(I)180:                   PIC(II)=PIC(I)
179:                ENDIF181:                ENDIF
180:             ENDIF182:             ENDIF
181:          ENDDO183:          ENDDO
182:          LENIC=II184:          LENIC=II
183:       ENDIF185:       ENDIF
184: C186: C
185:       IF(LFILL) THEN187:       IF(LFILL) THEN
 188:          PRINT*, "msb50 in l1081 intcor2.f"
186:          PRINT*, "LENIC", LENIC189:          PRINT*, "LENIC", LENIC
187:          DO II=1,LENIC190:          DO II=1,LENIC
188:             IF(B1IC(II).LE.0.001 .OR.LALL) THEN191:             IF(B1IC(II).LE.0.001 .OR.LALL) THEN
189:                II1=KAR(II)192:                II1=KAR(II)
190:                JJ1=LAR(II)193:                JJ1=LAR(II)
 194:                PRINT*, "msb50, B1IC, II1, JJ1", II1, JJ1
191:                IF(II1.GT.0.AND.JJ1.GT.0) THEN195:                IF(II1.GT.0.AND.JJ1.GT.0) THEN
192:                   QAT(1)=II1196:                   QAT(1)=II1
193:                   QAT(2)=JJ1197:                   QAT(2)=JJ1
194:                   CALL CODES(QAT(3),0,0,0,198:                   CALL CODES(QAT(3),0,0,0,
195:      &                       NATOM,IMOVE,IAC,1,QAT(1),QAT(2),199:      &                       NATOM,IMOVE,IAC,1,QAT(1),QAT(2),
196:      &                       0,0,0,0,0,0,0,0,0,0,0,0,0,0,200:      &                       0,0,0,0,0,0,0,0,0,0,0,0,0,0,
197:      &                       .FALSE.,0,                   ! DRUDE201:      &                       .FALSE.,0,                   ! DRUDE
198:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP202:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP
199:      &                       .TRUE.,.TRUE.)203:      &                       .TRUE.,.TRUE.)
200:                   IX=QAT(3)204:                   IX=QAT(3)
 205:                   PRINT*, "msb50 B1IC IX", IX
201:                   IF(IX.GT.0) THEN206:                   IF(IX.GT.0) THEN
202:                      B1IC(II)=CBB(IX)207:                      B1IC(II)=CBB(IX)
 208:                      PRINT*, "msb50 B1IC II",II, B1IC(II)
203:                   ELSE209:                   ELSE
204:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(5A)')210:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(5A)')
205:      &               ' %BUILDC-WARNING: no bond parameters for types (',211:      &               ' %BUILDC-WARNING: no bond parameters for types (',
206:      &                ATC(IAC(II1)),',',ATC(IAC(JJ1)),')'212:      &                ATC(IAC(II1)),',',ATC(IAC(JJ1)),')'
207:                   ENDIF213:                   ENDIF
208:                ENDIF214:                ENDIF
209:             ENDIF215:             ENDIF
210: C        216: C        
211:             IF(B2IC(II).LE.0.001 .OR.LALL) THEN217:             IF(B2IC(II).LE.0.001 .OR.LALL) THEN
212:                II1=IAR(II)218:                II1=IAR(II)
215:                IF(II1.GT.0.AND.JJ1.GT.0) THEN221:                IF(II1.GT.0.AND.JJ1.GT.0) THEN
216:                   QAT(1)=II1222:                   QAT(1)=II1
217:                   QAT(2)=JJ1223:                   QAT(2)=JJ1
218:                   CALL CODES(QAT(3),0,0,0,224:                   CALL CODES(QAT(3),0,0,0,
219:      &                       NATOM,IMOVE,IAC,1,QAT(1),QAT(2),225:      &                       NATOM,IMOVE,IAC,1,QAT(1),QAT(2),
220:      &                       0,0,0,0,0,0,0,0,0,0,0,0,0,0,226:      &                       0,0,0,0,0,0,0,0,0,0,0,0,0,0,
221:      &                       .FALSE.,0,                   ! DRUDE227:      &                       .FALSE.,0,                   ! DRUDE
222:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP228:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP
223:      &                       .TRUE.,.TRUE.)229:      &                       .TRUE.,.TRUE.)
224:                   IX=QAT(3)230:                   IX=QAT(3)
 231:                   PRINT*, "msb50, B2IC, IX", IX
225:                   IF(IX.GT.0) THEN232:                   IF(IX.GT.0) THEN
226:                      B2IC(II)=CBB(IX)233:                      B2IC(II)=CBB(IX)
 234:                      PRINT*, "msb50 B2IC II", II, B2IC(II)
227:                   ELSE235:                   ELSE
228:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(5A)')236:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(5A)')
229:      &               ' %BUILDC-WARNING: no bond parameters for types (',237:      &               ' %BUILDC-WARNING: no bond parameters for types (',
230:      &                ATC(IAC(II1)),',',ATC(IAC(JJ1)),')'238:      &                ATC(IAC(II1)),',',ATC(IAC(JJ1)),')'
231:                   ENDIF239:                   ENDIF
232:                ENDIF240:                ENDIF
233:             ENDIF241:             ENDIF
234: C        242: C        
235:             IF(T1IC(II).LE.0.001 .OR.LALL) THEN243:             IF(T1IC(II).LE.0.001 .OR.LALL) THEN
236:                II1=LAR(II)244:                II1=LAR(II)
241:                   QAT(2)=JJ1249:                   QAT(2)=JJ1
242:                   QAT(3)=KK1250:                   QAT(3)=KK1
243:                   CALL CODES(0,QAT(4),0,0,251:                   CALL CODES(0,QAT(4),0,0,
244:      &                       NATOM,IMOVE,IAC,0,0,0,252:      &                       NATOM,IMOVE,IAC,0,0,0,
245:      &                       1,QAT(1),QAT(2),QAT(3),253:      &                       1,QAT(1),QAT(2),QAT(3),
246:      &                       0,0,0,0,0,0,0,0,0,0,254:      &                       0,0,0,0,0,0,0,0,0,0,
247:      &                       .FALSE.,0,                   ! DRUDE255:      &                       .FALSE.,0,                   ! DRUDE
248:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP256:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP
249:      &                       .TRUE.,.TRUE.)257:      &                       .TRUE.,.TRUE.)
250:                   IX=QAT(4)258:                   IX=QAT(4)
 259:                   PRINT*, "msb50 T1IC,",II1,JJ1,KK1, IX
251:                   IF(IX.GT.0) THEN260:                   IF(IX.GT.0) THEN
252:                      T1IC(II)=CTB(IX)*RADDEG261:                      T1IC(II)=CTB(IX)*RADDEG
 262:                      PRINT*, "msb50 T1IC", T1IC(II)
253:                   ELSE263:                   ELSE
254:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(7A)')264:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(7A)')
255:      &              ' %BUILDC-WARNING: no angle parameters for types (',265:      &              ' %BUILDC-WARNING: no angle parameters for types (',
256:      &             ATC(IAC(II1)),',',ATC(IAC(JJ1)),',',ATC(IAC(KK1)),')'266:      &             ATC(IAC(II1)),',',ATC(IAC(JJ1)),',',ATC(IAC(KK1)),')'
257:                   ENDIF267:                   ENDIF
258:                ENDIF268:                ENDIF
259:             ENDIF269:             ENDIF
260: C        270: C        
261:             IF(T2IC(II).LE.0.001 .OR.LALL) THEN271:             IF(T2IC(II).LE.0.001 .OR.LALL) THEN
262:                II1=IAR(II)272:                II1=IAR(II)
271:                   QAT(2)=JJ1281:                   QAT(2)=JJ1
272:                   QAT(3)=KK1282:                   QAT(3)=KK1
273:                   CALL CODES(0,QAT(4),0,0,283:                   CALL CODES(0,QAT(4),0,0,
274:      &                       NATOM,IMOVE,IAC,0,0,0,284:      &                       NATOM,IMOVE,IAC,0,0,0,
275:      &                       1,QAT(1),QAT(2),QAT(3),285:      &                       1,QAT(1),QAT(2),QAT(3),
276:      &                       0,0,0,0,0,0,0,0,0,0,286:      &                       0,0,0,0,0,0,0,0,0,0,
277:      &                       .FALSE.,0,                   ! DRUDE287:      &                       .FALSE.,0,                   ! DRUDE
278:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP288:      &                       0,0,0,0,0,0,0,0,0,0,         !##CMAP
279:      &                       .TRUE.,.TRUE.)289:      &                       .TRUE.,.TRUE.)
280:                   IX=QAT(4)290:                   IX=QAT(4)
 291:                   !PRINT*, "msb50 T2IC", II1,JJ1,KK1,IX
281:                   IF(IX.GT.0) THEN292:                   IF(IX.GT.0) THEN
282:                      T2IC(II)=CTB(IX)*RADDEG293:                      T2IC(II)=CTB(IX)*RADDEG
 294:                      PRINT*, "msb50 T2IC",T2IC(II) 
283:                   ELSE295:                   ELSE
284:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(7A)')296:                      IF(WRNLEV.GE.2) WRITE(OUTU,'(7A)')
285:      &              ' %BUILDC-WARNING: no angle parameters for types (',297:      &              ' %BUILDC-WARNING: no angle parameters for types (',
286:      &             ATC(IAC(II1)),',',ATC(IAC(JJ1)),',',ATC(IAC(KK1)),')'298:      &             ATC(IAC(II1)),',',ATC(IAC(JJ1)),',',ATC(IAC(KK1)),')'
287:                   ENDIF299:                   ENDIF
288:                ENDIF300:                ENDIF
289:             ENDIF301:             ENDIF
290: C302: C
291: C           fill dihedrals (if explicitly requested)303: C           fill dihedrals (if explicitly requested)
292:             IF(LDIHE .AND. .NOT.TAR(II)) THEN304:             IF(LDIHE .AND. .NOT.TAR(II)) THEN
 305:               PRINT*, "msb50, in dihedral run"
293:               IF(PIC(II).LE.0.001 .OR.LALL) THEN306:               IF(PIC(II).LE.0.001 .OR.LALL) THEN
294:                 II1=IAR(II)307:                 II1=IAR(II)
295:                 JJ1=JAR(II)308:                 JJ1=JAR(II)
296:                 KK1=KAR(II)309:                 KK1=KAR(II)
297:                 LL1=LAR(II)310:                 LL1=LAR(II)
298:                 IF(II1.GT.0.AND.JJ1.GT.0.AND.KK1.GT.0.AND.LL1.GT.0) THEN311:                 IF(II1.GT.0.AND.JJ1.GT.0.AND.KK1.GT.0.AND.LL1.GT.0) THEN
299:                   QAT(1)=II1312:                   QAT(1)=II1
300:                   QAT(2)=JJ1313:                   QAT(2)=JJ1
301:                   QAT(3)=KK1314:                   QAT(3)=KK1
302:                   QAT(4)=LL1315:                   QAT(4)=LL1


r15547/intcor.src 2017-01-21 10:33:07.847188235 +0000 r15546/intcor.src 2017-01-21 10:33:15.174250000 +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"
100:       IF(WRD.EQ.'READ') THEN101:       IF(WRD.EQ.'READ') THEN
101: C  READ-INTERNAL-COORDINATE-FILE102: C  READ-INTERNAL-COORDINATE-FILE
102:          IUNIC=GTRMI(COMLYN,COMLEN,'UNIT',-1)103:          IUNIC=GTRMI(COMLYN,COMLEN,'UNIT',-1)
103:          LAPPE=INDXA(COMLYN,COMLEN,'APPE').NE.0104:          LAPPE=INDXA(COMLYN,COMLEN,'APPE').NE.0
104:          ICARD=1105:          ICARD=1
105:          IF(INDXA(COMLYN,COMLEN,'FILE').NE.0) ICARD=0106:          IF(INDXA(COMLYN,COMLEN,'FILE').NE.0) ICARD=0
106:          IF(INDXA(COMLYN,COMLEN,'CARD').NE.0) ICARD=1107:          IF(INDXA(COMLYN,COMLEN,'CARD').NE.0) ICARD=1
107:          ISTART=1108:          ISTART=1
108:          IF(ICARD.EQ.1 .AND. IUNIC.LT.0) IUNIC=ISTRM109:          IF(ICARD.EQ.1 .AND. IUNIC.LT.0) IUNIC=ISTRM
109:          IF(INDXA(COMLYN,COMLEN,'SAVE').NE.0) THEN110:          IF(INDXA(COMLYN,COMLEN,'SAVE').NE.0) THEN
402:              I=LENIC+100  ! free excess space403:              I=LENIC+100  ! free excess space
403:              CALL REINTC(I,INTLEN,LENIC,BINTCR,LINTCR)404:              CALL REINTC(I,INTLEN,LENIC,BINTCR,LINTCR)
404:            ENDIF405:            ENDIF
405:          ENDIF406:          ENDIF
406: C-----------------------------------------------------------------------407: C-----------------------------------------------------------------------
407:       ELSE IF(WRD.EQ.'PARA') THEN408:       ELSE IF(WRD.EQ.'PARA') THEN
408: C  PROCESS-PARAMETER-COMMAND409: C  PROCESS-PARAMETER-COMMAND
409:          LAPPE=INDXA(COMLYN,COMLEN,'ALL').NE.0410:          LAPPE=INDXA(COMLYN,COMLEN,'ALL').NE.0
410:          LDIHE=INDXA(COMLYN,COMLEN,'DIHE').NE.0411:          LDIHE=INDXA(COMLYN,COMLEN,'DIHE').NE.0
411:          LIMPR=INDXA(COMLYN,COMLEN,'IMPR').NE.0412:          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
412:          CALL PURGIC(LENIC,417:          CALL PURGIC(LENIC,
413:      &               HEAP(BINTCR(INTB1)),HEAP(BINTCR(INTB2)),418:      &               HEAP(BINTCR(INTB1)),HEAP(BINTCR(INTB2)),
414:      &               HEAP(BINTCR(INTT1)),HEAP(BINTCR(INTT2)),419:      &               HEAP(BINTCR(INTT1)),HEAP(BINTCR(INTT2)),
415:      &               HEAP(BINTCR(INTPIC)),HEAP(BINTCR(INTIAR)),420:      &               HEAP(BINTCR(INTPIC)),HEAP(BINTCR(INTIAR)),
416:      &               HEAP(BINTCR(INTJAR)),HEAP(BINTCR(INTKAR)),421:      &               HEAP(BINTCR(INTJAR)),HEAP(BINTCR(INTKAR)),
417:      &               HEAP(BINTCR(INTLAR)),HEAP(BINTCR(INTTAR)),422:      &               HEAP(BINTCR(INTLAR)),HEAP(BINTCR(INTTAR)),
418:      &               .TRUE.,LAPPE,LDIHE,LIMPR,.FALSE.,423:      &               .TRUE.,LAPPE,LDIHE,LIMPR,.FALSE.,
419:      &               NATOM,IMOVE,IAC)424:      &               NATOM,IMOVE,IAC)
420: C-----------------------------------------------------------------------425: C-----------------------------------------------------------------------
421:       ELSE IF(WRD.EQ.'DIFF') THEN426:       ELSE IF(WRD.EQ.'DIFF') THEN


r15547/keywords.f 2017-01-21 10:33:12.814250000 +0000 r15546/keywords.f 2017-01-21 10:33:19.538250000 +0000
 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,
 29:      & nocistransrna,nocistransdna,atmass1,checkcistransalways,AMBERICT,AMBSTEPT,AMBIT, 29:      & nocistransrna,nocistransdna,atmass1,checkcistransalways,AMBERICT,AMBSTEPT,AMBIT,
 30:      & AMBPERTT, PERTHRESH, AMBOLDPERTT,AMBICDNEBT,  30:      & AMBPERTT, PERTHRESH, AMBOLDPERTT, 
 31:      & ambpdb_unit, ambrst_unit, myunitnew, mdcrd_unit, mdinfo_unit 31:      & ambpdb_unit, ambrst_unit, myunitnew, mdcrd_unit, mdinfo_unit
 32:       USE MODNEB 32:       USE MODNEB
 33:       USE modmxatms   ! needed for CHARMM 33:       USE modmxatms   ! needed for CHARMM
 34:       USE modcharmm    34:       USE modcharmm   
 35:       USE MODUNRES 35:       USE MODUNRES
 36:       USE KEYNEB, NNNIMAGE=>NIMAGE 36:       USE KEYNEB, NNNIMAGE=>NIMAGE
 37:       USE KEYCONNECT 37:       USE KEYCONNECT
 38:       USE MODMEC 38:       USE MODMEC
 39:       USE MODGUESS 39:       USE MODGUESS
 40:       USE PORFUNCS 40:       USE PORFUNCS
539:       BISECTMINDIST=1.0539:       BISECTMINDIST=1.0
540:       BISECTMAXATTEMPTS=5540:       BISECTMAXATTEMPTS=5
541:       CHRIGIDT=.FALSE.541:       CHRIGIDT=.FALSE.
542:       PTRANS=0.D0542:       PTRANS=0.D0
543:       TRANSMAX=0.D0543:       TRANSMAX=0.D0
544:       PROT=0.D0544:       PROT=0.D0
545:       ROTMAX=0.D0545:       ROTMAX=0.D0
546:       BBRSDMT=.FALSE.546:       BBRSDMT=.FALSE.
547:       AMBERICT=.FALSE.547:       AMBERICT=.FALSE.
548:       AMBSTEPT=.FALSE.548:       AMBSTEPT=.FALSE.
549:       AMBICDNEBT = .FALSE. 
550:       AMBPERTT=.FALSE.549:       AMBPERTT=.FALSE.
551:       AMBOLDPERTT=.FALSE.550:       AMBOLDPERTT=.FALSE.
552:       AMBIT = .FALSE. 551:       AMBIT = .FALSE. 
553: 552: 
554:       DIJKSTRALOCAL=1.0D0553:       DIJKSTRALOCAL=1.0D0
555:       DNEBSWITCH=-1.0D0554:       DNEBSWITCH=-1.0D0
556:       PATHSDSTEPS=-1555:       PATHSDSTEPS=-1
557:       NUSEEV=-1556:       NUSEEV=-1
558: ! sf344557: ! sf344
559:       PYGPERIODICT=.FALSE.558:       PYGPERIODICT=.FALSE.
748:  747:  
749:       ELSE IF (WORD.eq.'AMBPERTOLD') THEN748:       ELSE IF (WORD.eq.'AMBPERTOLD') THEN
750:         PRINT*, "original perturbation scheme"749:         PRINT*, "original perturbation scheme"
751:         AMBOLDPERTT = .TRUE.750:         AMBOLDPERTT = .TRUE.
752: 751: 
753:       ELSE IF (WORD.eq. 'AMBPERTONLY') THEN752:       ELSE IF (WORD.eq. 'AMBPERTONLY') THEN
754:         AMBPERTT = .TRUE.753:         AMBPERTT = .TRUE.
755:         CALL READF(PERTHRESH)        754:         CALL READF(PERTHRESH)        
756:         PRINT*, "amber pertonly, perthresh", perthresh755:         PRINT*, "amber pertonly, perthresh", perthresh
757: 756: 
758:       ELSE IF (WORD.eq. 'AMBICDNEB') THEN 
759:         AMBICDNEBT = .TRUE. 
760:  
761:       ELSE IF (WORD.eq.'NAB') THEN757:       ELSE IF (WORD.eq.'NAB') THEN
762:         IF (FREEZERES) CALL A9RESTOATOM(FROZENRES,FROZEN,NFREEZE)758:         IF (FREEZERES) CALL A9RESTOATOM(FROZENRES,FROZEN,NFREEZE)
763:         NABT=.TRUE.759:         NABT=.TRUE.
764:         DO i=1,3*NATOMS760:         DO i=1,3*NATOMS
765:                 Q(i) = COORDS1(i)761:                 Q(i) = COORDS1(i)
766:         END DO762:         END DO
767: ! save atom names in array zsym763: ! save atom names in array zsym
768:         do i=1,natoms764:         do i=1,natoms
769:                 zsym(i) = ih(m04+i-1)765:                 zsym(i) = ih(m04+i-1)
770:         end do766:         end do
2618:                ENDIF2614:                ENDIF
2619:                IF (NSWAP(J1).GT.3) THEN2615:                IF (NSWAP(J1).GT.3) THEN
2620:                   PRINT '(2(A,I8))','keyword> ERROR - number of swaps ',NSWAP(J1),' is > 3'2616:                   PRINT '(2(A,I8))','keyword> ERROR - number of swaps ',NSWAP(J1),' is > 3'
2621:                   STOP2617:                   STOP
2622:                ENDIF2618:                ENDIF
2623:                IF (NDUMMY+NPERMSIZE(J1)-1.GT.NATOMS) THEN2619:                IF (NDUMMY+NPERMSIZE(J1)-1.GT.NATOMS) THEN
2624:                   PRINT '(2(A,I8))','keyword> ERROR - number of atoms to be permuted in all groups is > number of atoms'2620:                   PRINT '(2(A,I8))','keyword> ERROR - number of atoms to be permuted in all groups is > number of atoms'
2625:                   STOP2621:                   STOP
2626:                ENDIF2622:                ENDIF
2627:                READ(1,*) PERMGROUP(NDUMMY:NDUMMY+NPERMSIZE(J1)-1),(SWAP1(J1,J2),SWAP2(J1,J2),J2=1,NSWAP(J1))2623:                READ(1,*) PERMGROUP(NDUMMY:NDUMMY+NPERMSIZE(J1)-1),(SWAP1(J1,J2),SWAP2(J1,J2),J2=1,NSWAP(J1))
2628:              !msb50 
2629:              !PRINT*,"PERMGROUP",PERMGROUP(NDUMMY:NDUMMY+NPERMSIZE(J1)-1),'swap',SWAP1(J1,1:NSWAP(J1)),SWAP2(J1,1:NSWAP(J1)) 
2630:                NDUMMY=NDUMMY+NPERMSIZE(J1)2624:                NDUMMY=NDUMMY+NPERMSIZE(J1)
2631:             ENDDO2625:             ENDDO
2632: !2626: !
2633: !  And another sanity check!2627: !  And another sanity check!
2634: !2628: !
2635:             DO J1=1,NDUMMY2629:             DO J1=1,NDUMMY
2636:                DO J2=J1+1,NDUMMY2630:                DO J2=J1+1,NDUMMY
2637:                   IF (PERMGROUP(J2).EQ.PERMGROUP(J1)) THEN2631:                   IF (PERMGROUP(J2).EQ.PERMGROUP(J1)) THEN
2638:                      PRINT '(2(A,I8))','keyword> ERROR - atom ',PERMGROUP(J1),' appears more than once'2632:                      PRINT '(2(A,I8))','keyword> ERROR - atom ',PERMGROUP(J1),' appears more than once'
2639:                      STOP2633:                      STOP


r15547/Makefile 2017-01-21 10:33:09.970250000 +0000 r15546/Makefile 2017-01-21 10:33:16.834250000 +0000
 11: #   but WITHOUT ANY WARRANTY; without even the implied warranty of 11: #   but WITHOUT ANY WARRANTY; without even the implied warranty of
 12: #   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: #   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: #   GNU General Public License for more details. 13: #   GNU General Public License for more details.
 14: # 14: #
 15: #   You should have received a copy of the GNU General Public License 15: #   You should have received a copy of the GNU General Public License
 16: #   along with this program; if not, write to the Free Software 16: #   along with this program; if not, write to the Free Software
 17: #   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: #   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: # 18: #
 19: OBJS = potential.o key.o commons.o modtwoend.o modhess.o modneb.o modamber.o modamber2.o modamber9.o\ 19: OBJS = potential.o key.o commons.o modtwoend.o modhess.o modneb.o modamber.o modamber2.o modamber9.o\
 20: modunres.o modguess.o modmec.o syminf.o vecck.o zwk.o \ 20: modunres.o modguess.o modmec.o syminf.o vecck.o zwk.o \
 21: getparams.o adm.o amb_natinterns.o capsid.o axdiff.o axpairs.o aziz.o c60diff.o c60p.o beig.o bfgsts.o \ 21: getparams.o adm.o capsid.o axdiff.o axpairs.o aziz.o c60diff.o c60p.o beig.o bfgsts.o \
 22: dimer.o guesspath.o unguesspath.o \ 22: dimer.o guesspath.o unguesspath.o \
 23: charmmBildc.o dtrap.o dumpit.o dumpp.o eig.o eigensort.o emie.o escp.o etrap.o fetchz.o geopt.o gmetry.o hessout.o \ 23: charmmBildc.o dtrap.o dumpit.o dumpp.o eig.o eigensort.o emie.o escp.o etrap.o fetchz.o geopt.o gmetry.o hessout.o \
 24: ions.o JM.o keywords.o h2o.o latmin.o ljdiff.o ljpdiff.o mdiff.o mied.o miel.o mlatmin.o \ 24: ions.o JM.o keywords.o h2o.o latmin.o ljdiff.o ljpdiff.o mdiff.o mied.o miel.o mlatmin.o \
 25: morse.o rotd.o mpdiff.o msdiff.o mslatmin.o oldneb.o pertable.o shifth.o dprand.o \ 25: morse.o rotd.o mpdiff.o msdiff.o mslatmin.o oldneb.o pertable.o shifth.o dprand.o \
 26: rotcon.o rotm.o scdiff.o scl.o secdiag.o siaz.o sortc.o sortxyz.o symmetry.o tcheck.o vdump.o \ 26: rotcon.o rotm.o scdiff.o scl.o secdiag.o siaz.o sortc.o sortxyz.o symmetry.o tcheck.o vdump.o \
 27: efol.o utils.o ljpbin.o ljpshift.o SW.o ctest.o double.o Clatmin.o input.o cpmdlatmin.o ljpkob.o \ 27: efol.o utils.o ljpbin.o ljpshift.o SW.o ctest.o double.o Clatmin.o input.o cpmdlatmin.o ljpkob.o \
 28: changep.o twoend.o dzugutov.o Zetterling.o rad.o odesd.o g2special.o PV.o PachecoC60.o xmylbfgs.o \ 28: changep.o twoend.o dzugutov.o Zetterling.o rad.o odesd.o g2special.o PV.o PachecoC60.o xmylbfgs.o \
 29: mylbfgs.o 2Dfunc.o OPTIM.o path.o connect.o mindist.o repelsp.o TIPmodes.o caardiff.o mycpu_time.o \ 29: mylbfgs.o 2Dfunc.o OPTIM.o path.o connect.o mindist.o repelsp.o TIPmodes.o caardiff.o mycpu_time.o \
 30: g46merdiff.o p46merdiff.o MB.o c10.o ectrap.o dctrap.o welch.o tosifumi.o tight.o dftb.o \ 30: g46merdiff.o p46merdiff.o MB.o c10.o ectrap.o dctrap.o welch.o tosifumi.o tight.o dftb.o \
 31: fd.o ljms.o meccano.o unmeccano.o porfuncs.o modcharmm.o modmxatms.o rigidfuncs.o tip.o \ 31: fd.o ljms.o meccano.o unmeccano.o porfuncs.o modcharmm.o modmxatms.o rigidfuncs.o tip.o \
 32: msevb_common.o msevb_interfaces.o checkedtrig.o msevb.o fmsevb.o drvmsevb.o cleanmemory.o Natb.o \ 32: msevb_common.o msevb_interfaces.o checkedtrig.o msevb.o fmsevb.o drvmsevb.o cleanmemory.o Natb.o \
 33: header.o wc.o binaryio.o BLN.o morph.o newmindist.o thomson.o minpermdist.o minperm.o inertia.o \ 33: header.o wc.o binaryio.o BLN.o morph.o newmindist.o thomson.o minpermdist.o minperm.o inertia.o \
 34: greatcirc.o GoOptim.o cubspl.o specialfuncts.o gsdata.o cubsplstring.o dqag.o \ 34: greatcirc.o GoOptim.o cubspl.o specialfuncts.o gsdata.o cubsplstring.o dqag.o \
 35: intcommons.o intcoords.o SiO2.o stock.o bhinterp.o Z2faster.o p4diff.o Gupta.o myorient.o hybridmin.o \ 35: intcommons.o intcoords.o SiO2.o stock.o bhinterp.o Z2faster.o p4diff.o Gupta.o myorient.o hybridmin.o \
 36: bbrsdm.o projh.o bisect.o dbpg.o lwotpgh.o gb.o gbbox.o gbd.o rigidb.o newcapsid.o newtip.o pyg.o \ 36: bbrsdm.o projh.o bisect.o dbpg.o lwotpgh.o gb.o gbbox.o gbd.o rigidb.o newcapsid.o newtip.o pyg.o stockghaa.o gay-berne.o growstring.o
 37: stockghaa.o gay-berne.o growstring.o ambnatintern_extras.o 
 38:  37: 
 39: GENF90FILES = porfuncs.f90 header.f90 38: GENF90FILES = porfuncs.f90 header.f90
 40:  39: 
 41: CHDUM = chdummy.o 40: CHDUM = chdummy.o
 42: UNDUM = unresdummy.o 41: UNDUM = unresdummy.o
 43: AMDUM = amberdummy.o 42: AMDUM = amberdummy.o
 44: AMHDUM = amhdummy.o 43: AMHDUM = amhdummy.o
 45:  44: 
 46: VPATH = NEB:CONNECT:CHARMM:AMH:. 45: VPATH = NEB:CONNECT:CHARMM:AMH:.
 47: LDFLAGS = -L. 46: LDFLAGS = -L.
136: ${unoptim}: $(AMB9DUM) $(AMHDUM) $(OBJS) $(EXTRAS) $(EXTRASUNRES) $(CHDUM) $(AMDUM)135: ${unoptim}: $(AMB9DUM) $(AMHDUM) $(OBJS) $(EXTRAS) $(EXTRASUNRES) $(CHDUM) $(AMDUM)
137:         $(FC) $(FFLAGS) ${SEARCH_PATH} -o ${unoptim} $(EXTRAS) $(EXTRASUNRES) $(CHDUM) $(AMDUM) $(AMHDUM) $(AMB9DUM) $(OBJS) \136:         $(FC) $(FFLAGS) ${SEARCH_PATH} -o ${unoptim} $(EXTRAS) $(EXTRASUNRES) $(CHDUM) $(AMDUM) $(AMHDUM) $(AMB9DUM) $(OBJS) \
138:         $(LDFLAGS) $(LIBS)137:         $(LDFLAGS) $(LIBS)
139:         rm libmyblas.a    libmylapack.a138:         rm libmyblas.a    libmylapack.a
140: 139: 
141: ${amoptim}: $(AMB9DUM) $(AMHDUM) $(OBJS) $(EXTRAS) $(EXTRASAMBER) $(CHDUM) $(UNDUM) 140: ${amoptim}: $(AMB9DUM) $(AMHDUM) $(OBJS) $(EXTRAS) $(EXTRASAMBER) $(CHDUM) $(UNDUM) 
142:         $(FC) $(FFLAGS) ${SEARCH_PATH} -o ${amoptim} $(EXTRAS) $(EXTRASAMBER) $(CHDUM) $(UNDUM) $(AMHDUM) $(AMB9DUM) $(OBJS) \141:         $(FC) $(FFLAGS) ${SEARCH_PATH} -o ${amoptim} $(EXTRAS) $(EXTRASAMBER) $(CHDUM) $(UNDUM) $(AMHDUM) $(AMB9DUM) $(OBJS) \
143:         $(LDFLAGS) $(LIBS)142:         $(LDFLAGS) $(LIBS)
144:         rm libmyblas.a    libmylapack.a143:         rm libmyblas.a    libmylapack.a
145: 144: 
146: ${amb9optim}: $(AMHDUM) $(EXTRAS) $(OBJS) $(CHDUM) $(AMDUM) $(UNDUM) libamber.a libnab.a 145: ${amb9optim}: $(AMHDUM) $(OBJS) $(CHDUM) $(AMDUM) $(UNDUM) libamber.a libnab.a 
147:         $(FC) $(FFLAGS) ${SEARCH_PATH} -o ${amb9optim} $(CHDUM) $(UNDUM) $(AMHDUM) $(AMDUM) $(EXTRAS) $(OBJS) \146:         $(FC) $(FFLAGS) ${SEARCH_PATH} -o ${amb9optim} $(CHDUM) $(UNDUM) $(AMHDUM) $(AMDUM) $(OBJS) \
148:         $(LDFLAGS) $(LIBS) libamber.a libnab.a $(LIBS)147:         $(LDFLAGS) $(LIBS) libamber.a libnab.a $(LIBS)
149: #        rm -f libmyblas.a    libmylapack.a modamber9.o modamber9.mod148: #        rm -f libmyblas.a    libmylapack.a modamber9.o modamber9.mod
150: 149: 
151: ${amhoptim}: $(AMB9DUM) $(EXTRASAMH) $(OBJS) $(EXTRAS) $(CHDUM) $(AMDUM) $(UNDUM) libamh.a 150: ${amhoptim}: $(AMB9DUM) $(EXTRASAMH) $(OBJS) $(EXTRAS) $(CHDUM) $(AMDUM) $(UNDUM) libamh.a 
152:         $(FC) $(FFLAGS) ${SEARCH_PATH} -o ${amhoptim} $(EXTRAS) $(EXTRASAMH)  $(CHDUM) $(UNDUM) $(AMDUM) $(AMB9DUM) $(OBJS) \151:         $(FC) $(FFLAGS) ${SEARCH_PATH} -o ${amhoptim} $(EXTRAS) $(EXTRASAMH)  $(CHDUM) $(UNDUM) $(AMDUM) $(AMB9DUM) $(OBJS) \
153:         $(LDFLAGS) $(LIBS)  $(LIBS) libamh.a152:         $(LDFLAGS) $(LIBS)  $(LIBS) libamh.a
154:         rm libmyblas.a    libmylapack.a 153:         rm libmyblas.a    libmylapack.a 
155: 154: 
156: timing:155: timing:
157:         rm -f $(OPROG) 156:         rm -f $(OPROG) 
204:         cd CHARMM; make FC="${FC}" FFLAGS="${FFLAGS} ${SEARCH_PATH}" PREFLX="${PREFLX}" PREFDIR="${PREFDIR}" \203:         cd CHARMM; make FC="${FC}" FFLAGS="${FFLAGS} ${SEARCH_PATH}" PREFLX="${PREFLX}" PREFDIR="${PREFDIR}" \
205:         CTYPE="${CTYPE}" FCMDIR=${FCMDIR} C31SRC="${C31SRC}" SRC31="${SRC31}"204:         CTYPE="${CTYPE}" FCMDIR=${FCMDIR} C31SRC="${C31SRC}" SRC31="${SRC31}"
206: SAT-Ghost:205: SAT-Ghost:
207: 206: 
208: ###################################### DEPENDENCIES ######################################207: ###################################### DEPENDENCIES ######################################
209: ${optim} ${coptim} ${unoptim} ${amoptim} ${amb9optim} ${amhoptim}: libnn.a libnc.a libmyblas.a libmylapack.a208: ${optim} ${coptim} ${unoptim} ${amoptim} ${amb9optim} ${amhoptim}: libnn.a libnc.a libmyblas.a libmylapack.a
210: ${amb9optim}: libamber.a modamber9.o209: ${amb9optim}: libamber.a modamber9.o
211: libamh.a: altpot_interfaces.o  amh_interfaces.o  globals_alt.o amhglobals.o210: libamh.a: altpot_interfaces.o  amh_interfaces.o  globals_alt.o amhglobals.o
212: libnn.a: modcharmm.o key.o commons.o modunres.o porfuncs.o modmec.o modguess.o efol.o growstring.o gsdata.o intcommons.o intcoords.o211: libnn.a: modcharmm.o key.o commons.o modunres.o porfuncs.o modmec.o modguess.o efol.o growstring.o gsdata.o intcommons.o intcoords.o
213: libnc.a: libnn.a key.o syminf.o modhess.o212: libnc.a: libnn.a key.o syminf.o modhess.o
214: libamber.a: commons.o modamber9.o 213: libamber.a: commons.o modamber9.o
215: libcharmm.a: key.o commons.o modtwoend.o modneb.o modhess.o modcharmm.o modmxatms.o intcommons.o214: libcharmm.a: key.o commons.o modtwoend.o modneb.o modhess.o modcharmm.o modmxatms.o intcommons.o
216: oldneb.o meccano.o unmeccano.o fetchz.o keywords.o OPTIM.o unresnebguessts.o: libnn.a215: oldneb.o meccano.o unmeccano.o fetchz.o keywords.o OPTIM.o unresnebguessts.o: libnn.a
217: fetchz.o keywords.o : libnc.a216: fetchz.o keywords.o : libnc.a
218: OPTIM.o connect.o: libnn.a libnc.a217: OPTIM.o connect.o: libnn.a libnc.a
219: 2Dfunc.o: modcharmm.o commons.o modhess.o218: 2Dfunc.o: modcharmm.o commons.o modhess.o
220: Clatmin.o: commons.o219: Clatmin.o: commons.o
221: JM.o: modhess.o porfuncs.o220: JM.o: modhess.o porfuncs.o
222: MB.o: commons.o modhess.o221: MB.o: commons.o modhess.o
223: OPTIM.o: key.o commons.o modhess.o modneb.o modtwoend.o vecck.o zwk.o modunres.o modmec.o modguess.o porfuncs.o 222: OPTIM.o: key.o commons.o modhess.o modneb.o modtwoend.o vecck.o zwk.o modunres.o modmec.o modguess.o porfuncs.o 
224: PV.o: commons.o key.o223: PV.o: commons.o key.o
225: PachecoC60.0: modhess.o224: PachecoC60.0: modhess.o
226: SW.o: key.o modhess.o225: SW.o: key.o modhess.o
227: TIPmodes.o: modhess.o226: TIPmodes.o: modhess.o
228: Zetterling.o: key.o modhess.o porfuncs.o227: Zetterling.o: key.o modhess.o porfuncs.o
229: adm.o: commons.o key.o porfuncs.o228: adm.o: commons.o key.o porfuncs.o
230: amber.o: commons.o key.o modamber.o modamber2.o modhess.o229: amber.o: commons.o key.o modamber.o modamber2.o modhess.o
231: amb_natinterns.o: commons.o intcommons.o modamber9.o 
232: ambnatintern_extras.o: commons.o intcommons.o modamber9.o myblas.o mylapack.o  
233: axdiff.o: modhess.o230: axdiff.o: modhess.o
234: aziz.o: modhess.o231: aziz.o: modhess.o
235: beig.o: commons.o key.o modtwoend.o modneb.o porfuncs.o232: beig.o: commons.o key.o modtwoend.o modneb.o porfuncs.o
236: bfgsts.o: commons.o key.o modtwoend.o vecck.o zwk.o modhess.o porfuncs.o233: bfgsts.o: commons.o key.o modtwoend.o vecck.o zwk.o modhess.o porfuncs.o
237: hybridmin.o: commons.o key.o modtwoend.o vecck.o zwk.o modhess.o porfuncs.o234: hybridmin.o: commons.o key.o modtwoend.o vecck.o zwk.o modhess.o porfuncs.o
238: c10.o: modhess.o235: c10.o: modhess.o
239: c60diff.o: modhess.o236: c60diff.o: modhess.o
240: c60p.o: modhess.o237: c60p.o: modhess.o
241: caardiff.o: modhess.o238: caardiff.o: modhess.o
242: capsid.o: modhess.o239: capsid.o: modhess.o


r15547/minpermdist.f90 2017-01-21 10:33:13.106250000 +0000 r15546/minpermdist.f90 2017-01-21 10:33:19.894250000 +0000
 64: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS) 64: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS)
 65: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), DX, DY, DZ, RMS, DBEST, XBEST(3*NATOMS) 65: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), DX, DY, DZ, RMS, DBEST, XBEST(3*NATOMS)
 66: DOUBLE PRECISION CMXA, CMXB, CMXC 66: DOUBLE PRECISION CMXA, CMXB, CMXC
 67: DOUBLE PRECISION ROTA(3,3), ROTINVA(3,3), ROTB(3,3), ROTINVB(3,3), ROTINVBBEST(3,3), ROTABEST(3,3), RMATBEST(3,3), TMAT(3,3) 67: DOUBLE PRECISION ROTA(3,3), ROTINVA(3,3), ROTB(3,3), ROTINVB(3,3), ROTINVBBEST(3,3), ROTABEST(3,3), RMATBEST(3,3), TMAT(3,3)
 68: DOUBLE PRECISION CMAX, CMAY, CMAZ, CMBX, CMBY, CMBZ, RMATCUMUL(3,3) 68: DOUBLE PRECISION CMAX, CMAY, CMAZ, CMBX, CMBY, CMBZ, RMATCUMUL(3,3)
 69: LOGICAL DEBUG, TWOD, RIGID, BULKT 69: LOGICAL DEBUG, TWOD, RIGID, BULKT
 70: DOUBLE PRECISION PDUMMYA(3*NATOMS), PDUMMYB(3*NATOMS), LDISTANCE, DUMMYC(3*NATOMS), XDUMMY 70: DOUBLE PRECISION PDUMMYA(3*NATOMS), PDUMMYB(3*NATOMS), LDISTANCE, DUMMYC(3*NATOMS), XDUMMY
 71:  71: 
 72: IF (INTMINPERMT.AND.(INTINTERPT.OR.DESMINT)) THEN 72: IF (INTMINPERMT.AND.(INTINTERPT.OR.DESMINT)) THEN
 73:    CALL INTMINPERM(COORDSB,COORDSA,DISTANCE,RMAT,DEBUG) 73:    CALL INTMINPERM(COORDSB,COORDSA,DISTANCE,RMAT,DEBUG)
 74:    PRINT*, "distance in minperm", distance 
 75:    RETURN 74:    RETURN
 76: ENDIF 75: ENDIF
 77:  76: 
 78: ! DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS) 77: ! DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS)
 79: ! DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS) 78: ! DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS)
 80:  79: 
 81: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 80: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 82: ! CALL OCHARMM(DUMMYA,VNEW,ENERGY,.FALSE.,.FALSE.) 81: ! CALL OCHARMM(DUMMYA,VNEW,ENERGY,.FALSE.,.FALSE.)
 83: ! CALL POTENTIAL(DUMMYA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 82: ! CALL POTENTIAL(DUMMYA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
 84: ! PRINT '(2(A,F25.15))',' Initial energy=',ENERGY,' RMS=',RMS 83: ! PRINT '(2(A,F25.15))',' Initial energy=',ENERGY,' RMS=',RMS
444:    COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!443:    COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!
445: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!444: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
446: !  DO J1=1,(NATOMS/2)445: !  DO J1=1,(NATOMS/2)
447: !     XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-COORDSA(3*(J1-1)+1))**2+ &446: !     XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-COORDSA(3*(J1-1)+1))**2+ &
448: ! &                 (COORDSB(3*(J1-1)+2)-COORDSA(3*(J1-1)+2))**2+ &447: ! &                 (COORDSB(3*(J1-1)+2)-COORDSA(3*(J1-1)+2))**2+ &
449: ! &                 (COORDSB(3*(J1-1)+3)-COORDSA(3*(J1-1)+3))**2448: ! &                 (COORDSB(3*(J1-1)+3)-COORDSA(3*(J1-1)+3))**2
450: !  ENDDO449: !  ENDDO
451: !  PRINT '(A,F20.10)','XDUMMY=',XDUMMY450: !  PRINT '(A,F20.10)','XDUMMY=',XDUMMY
452: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!451: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
453:       CLOSE(10)452:       CLOSE(10)
454:  
455:       IF (AMBERT) THEN 
456:          CALL check_valleu_chirality(COORDSB, COORDSA,DEBUG) 
457:          print*, "in minperm after check_valleu" 
458:       ENDIF 
459: RETURN453: RETURN
460: END SUBROUTINE MINPERMDIST454: END SUBROUTINE MINPERMDIST


r15547/modamber9.f90 2017-01-21 10:33:13.430250000 +0000 r15546/modamber9.f90 2017-01-21 10:33:20.194250000 +0000
 26:       CHARACTER(len=81) :: prmtop 26:       CHARACTER(len=81) :: prmtop
 27:       DOUBLE PRECISION, PARAMETER :: TEMPPARAM = 0.0000001 27:       DOUBLE PRECISION, PARAMETER :: TEMPPARAM = 0.0000001
 28: !      INTEGER   NATOM, irest, ntb 28: !      INTEGER   NATOM, irest, ntb
 29: !      DOUBLE PRECISION    tt,crd(3*656),vel(3,656) 29: !      DOUBLE PRECISION    tt,crd(3*656),vel(3,656)
 30:        INTEGER MYUNITNEW, MDINFO_UNIT, MDCRD_UNIT, AMBRST_UNIT, AMBPDB_UNIT 30:        INTEGER MYUNITNEW, MDINFO_UNIT, MDCRD_UNIT, AMBRST_UNIT, AMBPDB_UNIT
 31:        31:       
 32: !      CHARACTER(20) AMBERPRMTOP  32: !      CHARACTER(20) AMBERPRMTOP 
 33: !      DOUBLE PRECISION xbar, ybar, zbar 33: !      DOUBLE PRECISION xbar, ybar, zbar
 34: !   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE        :: y 34: !   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE        :: y
 35:    double precision,  dimension(:), allocatable :: x 35:    double precision,  dimension(:), allocatable :: x
 36:    integer, dimension(:), allocatable :: ipairs 36:    integer, dimension(:), allocatable :: ix, ipairs
 37:    integer, dimension(:), allocatable, target :: ix 
 38:    character(len=4), dimension(:), allocatable :: ih 37:    character(len=4), dimension(:), allocatable :: ih
 39:  38: 
 40: logical master 39: logical master
 41: common/extra_logical/master 40: common/extra_logical/master
 42:  41: 
 43:  42: 
 44: integer iredir(8) 43: integer iredir(8)
 45: common/nmrrdr/redir,iredir 44: common/nmrrdr/redir,iredir
 46:  45: 
 47: integer nmropt,iprint,noeskp,iscale,ipnlty,iuse,maxsub,jar 46: integer nmropt,iprint,noeskp,iscale,ipnlty,iuse,maxsub,jar
299: integer, dimension(:), allocatable :: NSIDECHAIN298: integer, dimension(:), allocatable :: NSIDECHAIN
300: integer, dimension(:), allocatable :: NCHIRAL299: integer, dimension(:), allocatable :: NCHIRAL
301: integer, dimension(:), allocatable :: PHIPSI300: integer, dimension(:), allocatable :: PHIPSI
302: integer, dimension(:), allocatable :: OMEGAC301: integer, dimension(:), allocatable :: OMEGAC
303: integer, dimension(:), allocatable :: TW_SIDECHAIN302: integer, dimension(:), allocatable :: TW_SIDECHAIN
304: integer, dimension(:), allocatable :: CHIRAL303: integer, dimension(:), allocatable :: CHIRAL
305: logical, dimension(:), allocatable :: IS_SIDECHAIN304: logical, dimension(:), allocatable :: IS_SIDECHAIN
306: integer, dimension(:), allocatable :: NICTOT !no of residues per segment305: integer, dimension(:), allocatable :: NICTOT !no of residues per segment
307: integer:: nseg306: integer:: nseg
308: logical:: tomegac307: logical:: tomegac
309: logical:: AMBERICT, AMBSTEPT, AMBPERTT,AMBIT, AMBOLDPERTT, AMBICDNEBT308: logical:: AMBERICT, AMBSTEPT, AMBPERTT,AMBIT, AMBOLDPERTT
310: double precision:: PERTHRESH !threshhold for dihedral perturbation309: double precision:: PERTHRESH !threshhold for dihedral perturbation
311:                     !perturb if difference (start-end) > perthresh310:                     !perturb if difference (start-end) > perthresh
312: integer:: NTW_ANGLES !if ambpertt: how many angles can you really twist311: integer:: NTW_ANGLES !if ambpertt: how many angles can you really twist
313: double precision, dimension(:), allocatable :: TW_DIFFP312: double precision, dimension(:), allocatable :: TW_DIFFP
314: !need this for natinterns313: 
315: integer :: NBOND314: 
316: logical, dimension (:),allocatable ::AM_BACKBONE 
317: integer, allocatable:: IBIB(:),JBJB(:)  
318: END MODULE modamber9 315: END MODULE modamber9 
319: 316: 


r15547/mylbfgs.f 2017-01-21 10:33:13.726250000 +0000 r15546/mylbfgs.f 2017-01-21 10:33:20.506250000 +0000
287:          ELSE IF (CHRMMT) THEN 287:          ELSE IF (CHRMMT) THEN 
288:             CALL GETKD(KD) ! get width of sparse band in G matrix KD288:             CALL GETKD(KD) ! get width of sparse band in G matrix KD
289:             CALL GETNNZ(NNZ) ! get number of non-zero elements in B-matrix289:             CALL GETNNZ(NNZ) ! get number of non-zero elements in B-matrix
290:             NOCOOR=.FALSE. ! calculate internals therefore NOCOOR is false290:             NOCOOR=.FALSE. ! calculate internals therefore NOCOOR is false
291:             NODERV = .FALSE.291:             NODERV = .FALSE.
292:             GINT(1:N)=0.0D0 ! to prevent NaN's for Sun!292:             GINT(1:N)=0.0D0 ! to prevent NaN's for Sun!
293:             XINT(1:N)=0.0D0 ! to prevent NaN's for Sun!293:             XINT(1:N)=0.0D0 ! to prevent NaN's for Sun!
294:             CALL TRANSFORM(X,G,XINT,GINT,N,3*NATOMS,NNZ,NOCOOR,NODERV,KD,INTEPSILON)294:             CALL TRANSFORM(X,G,XINT,GINT,N,3*NATOMS,NNZ,NOCOOR,NODERV,KD,INTEPSILON)
295:             OLDQ(1:N)=XINT(1:N)    ! store internals295:             OLDQ(1:N)=XINT(1:N)    ! store internals
296:             OLDGINT(1:N)=GINT(1:N) ! store gradient in internals296:             OLDGINT(1:N)=GINT(1:N) ! store gradient in internals
297:          ELSE IF (AMBERT) THEN 
298:             PRINT*, "not yet set up for amber" 
299:             STOP 
300:          ENDIF297:          ENDIF
301:       ENDIF298:       ENDIF
302: C299: C
303: C  for CHRMMT:300: C  for CHRMMT:
304: C  X       contains current Cartesians301: C  X       contains current Cartesians
305: C  G       contains current gradient302: C  G       contains current gradient
306: C  XINT    contains current internals303: C  XINT    contains current internals
307: C  GINT    contains current gradient in internals304: C  GINT    contains current gradient in internals
308: C  OLDQ    contains internals for initial geometry305: C  OLDQ    contains internals for initial geometry
309: C  OLDGINT contains gradient in internals for initial geometry306: C  OLDGINT contains gradient in internals for initial geometry


r15547/newneb.f90 2017-01-21 10:33:10.258250000 +0000 r15546/newneb.f90 2017-01-21 10:33:17.122250000 +0000
142:           X         => XYZ(NOPT+1:NOPT*(NIMAGE+1))142:           X         => XYZ(NOPT+1:NOPT*(NIMAGE+1))
143:           EIMAGE    => EEE(2:NIMAGE+1)143:           EIMAGE    => EEE(2:NIMAGE+1)
144:           G         => GGG(NOPT+1:NOPT*(NIMAGE+1))144:           G         => GGG(NOPT+1:NOPT*(NIMAGE+1))
145:           GSPR      => SSS(NOPT+1:NOPT*(NIMAGE+1))145:           GSPR      => SSS(NOPT+1:NOPT*(NIMAGE+1))
146:           RMSFIMAGE => RRR(2:NIMAGE+1)146:           RMSFIMAGE => RRR(2:NIMAGE+1)
147:           TANPTR => TANVEC147:           TANPTR => TANVEC
148:           EEE = 0.0D0148:           EEE = 0.0D0
149:           EEE(1)=EINITIAL149:           EEE(1)=EINITIAL
150:           EEE(NIMAGE+2)=EFINAL150:           EEE(NIMAGE+2)=EFINAL
151:           IF (DESMINT.AND..NOT.GROWSTRINGT) THEN151:           IF (DESMINT.AND..NOT.GROWSTRINGT) THEN
152:           PRINT*, "newneb>" 
153:              XYZCART(:3*NATOMS) = QQ152:              XYZCART(:3*NATOMS) = QQ
154:              XYZCART(3*NATOMS*(NIMAGE+1)+1:) = FINFIN153:              XYZCART(3*NATOMS*(NIMAGE+1)+1:) = FINFIN
155: 154: 
156:              PREVDIH => DIHINFO(1,:)155:              PREVDIH => DIHINFO(1,:)
157:              CALL CART2INT(QQ,XYZ(:NOPT))156:              CALL CART2INT(QQ,XYZ(:NOPT))
158: 157: 
159:              DO J1 = 2,NIMAGE+2158:              DO J1 = 2,NIMAGE+2
160:                 ! align all other dihedrals to start159:                 ! align all other dihedrals to start
161:                 DIHINFO(J1,:) = DIHINFO(1,:)160:                 DIHINFO(J1,:) = DIHINFO(1,:)
162:              ENDDO161:              ENDDO


r15547/nnutils.f90 2017-01-21 10:33:10.566250000 +0000 r15546/nnutils.f90 2017-01-21 10:33:17.402250000 +0000
 95:           DVEC=SQRT(DVEC) 95:           DVEC=SQRT(DVEC)
 96:           SEPARATION=SUM(DVEC) 96:           SEPARATION=SUM(DVEC)
 97:  97: 
 98:           DEVIATION = ABS( 100*( (NIMAGE+1)*DVEC/SEPARATION -1 ) ) 98:           DEVIATION = ABS( 100*( (NIMAGE+1)*DVEC/SEPARATION -1 ) )
 99:           AVDEV     = SUM(DEVIATION)/(NIMAGE+1) 99:           AVDEV     = SUM(DEVIATION)/(NIMAGE+1)
100:      END SUBROUTINE INTERNALIMAGEDISTRIBUTION100:      END SUBROUTINE INTERNALIMAGEDISTRIBUTION
101: 101: 
102:      SUBROUTINE MAKEIMAGE(EINITIAL,EFINAL,QQ,FINFIN) ! INTERPOLATES THE BAND102:      SUBROUTINE MAKEIMAGE(EINITIAL,EFINAL,QQ,FINFIN) ! INTERPOLATES THE BAND
103:           USE NEBDATA103:           USE NEBDATA
104:           USE KEYNEB,ONLY: NIMAGE104:           USE KEYNEB,ONLY: NIMAGE
105:           USE KEY,ONLY: MORPHT, MSTEPS, DEBUG, GREATCIRCLET, GCIMAGE, GCCONV, GCUPDATE, GCSTEPS, INTEPSILON, BULKT,AMBERT105:           USE KEY,ONLY: MORPHT, MSTEPS, DEBUG, GREATCIRCLET, GCIMAGE, GCCONV, GCUPDATE, GCSTEPS, INTEPSILON, BULKT
106:           USE INTCOMMONS, ONLY : NNZ, KD, NINTC, DESMINT, DIHINFO, PREVDIH, BACKTCUTOFF, INTERPBACKTCUT, MINBACKTCUT, &106:           USE INTCOMMONS, ONLY : NNZ, KD, NINTC, DESMINT, DIHINFO, PREVDIH, BACKTCUTOFF, INTERPBACKTCUT, MINBACKTCUT, &
107:                                  CHICDNEB107:                                  CHICDNEB
108:           USE COMMONS, ONLY : ZSYM, PARAM1,PARAM2,PARAM3108:           USE COMMONS, ONLY : ZSYM, PARAM1,PARAM2,PARAM3
109:           USE MODCHARMM, ONLY : CHRMMT, ICINTERPT109:           USE MODCHARMM, ONLY : CHRMMT, ICINTERPT
110:           USE MODAMBER9, ONLY: AMBERICT, AMBICDNEBT, NICTOT 
111:           USE PORFUNCS 110:           USE PORFUNCS 
112: 111: 
113:           IMPLICIT NONE112:           IMPLICIT NONE
114:           INTEGER :: I, J1, ITDONE, INTERVAL, NDONE, J2, K, ISTAT113:           INTEGER :: I, J1, ITDONE, INTERVAL, NDONE, J2, K, ISTAT
115:           DOUBLE PRECISION,ALLOCATABLE :: DELTAX(:)114:           DOUBLE PRECISION,ALLOCATABLE :: DELTAX(:)
116:           DOUBLE PRECISION DPRAND, SHIFT(NOPT), DUMMY, DLENGTH, DUMMY2, ENERGY, VNEW(NOPT), LRMS, &115:           DOUBLE PRECISION DPRAND, SHIFT(NOPT), DUMMY, DLENGTH, DUMMY2, ENERGY, VNEW(NOPT), LRMS, &
117:    &                       QS(NOPT), QF(NOPT), EINITIAL, EFINAL, COORDS(NOPT), SFRAC116:    &                       QS(NOPT), QF(NOPT), EINITIAL, EFINAL, COORDS(NOPT), SFRAC
118:           DOUBLE PRECISION, ALLOCATABLE :: XINITIAL(:), XIMAGE(:,:)117:           DOUBLE PRECISION, ALLOCATABLE :: XINITIAL(:), XIMAGE(:,:)
119:                     DOUBLE PRECISION,DIMENSION(:)         :: QQ,FINFIN118:                     DOUBLE PRECISION,DIMENSION(:)         :: QQ,FINFIN
120:           LOGICAL KNOWE, KNOWG, KNOWH, MFLAG119:           LOGICAL KNOWE, KNOWG, KNOWH, MFLAG
135: 134: 
136:           IF (MORPHT) THEN135:           IF (MORPHT) THEN
137:              IF (DESMINT) THEN136:              IF (DESMINT) THEN
138:                 print*, 'ERROR! MORPHT not implemented with DESM.'137:                 print*, 'ERROR! MORPHT not implemented with DESM.'
139:                 STOP138:                 STOP
140:              ENDIF139:              ENDIF
141: 140: 
142:              QS(1:3*NATOMS)=XYZ(1:NOPT)141:              QS(1:3*NATOMS)=XYZ(1:NOPT)
143:              QF(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)142:              QF(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)
144:              KNOWE=.FALSE.;  KNOWG=.FALSE.; KNOWH=.FALSE.143:              KNOWE=.FALSE.;  KNOWG=.FALSE.; KNOWH=.FALSE.
145:            
146:              CALL MORPH(MSTEPS,QS,QF,ENERGY,VNEW,MFLAG,LRMS,ITDONE,.TRUE.)144:              CALL MORPH(MSTEPS,QS,QF,ENERGY,VNEW,MFLAG,LRMS,ITDONE,.TRUE.)
147:              IF (ITDONE.LT.NIMAGE) THEN145:              IF (ITDONE.LT.NIMAGE) THEN
148:                 IF (MFLAG) THEN146:                 IF (MFLAG) THEN
149: !147: !
150: !  Unfortunately many statements in the NEB148: !  Unfortunately many statements in the NEB
151: !  routines do not specify the array bounds on array operations, which then default to the149: !  routines do not specify the array bounds on array operations, which then default to the
152: !  declared array size and therefore go wrong. Hence if we reduce the number of images here150: !  declared array size and therefore go wrong. Hence if we reduce the number of images here
153: !  then we have to deallocate and reallocate:151: !  then we have to deallocate and reallocate:
154: !                  152: !                  
155:                    PRINT '(A,I6)',' makeimage> Fewer morph points than images - reducing images to ',ITDONE153:                    PRINT '(A,I6)',' makeimage> Fewer morph points than images - reducing images to ',ITDONE
238:              ELSE236:              ELSE
239:                 DELTAX(1:NOPT) = ( XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2) ) - XYZ(1:NOPT) )/(NIMAGE+1)237:                 DELTAX(1:NOPT) = ( XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2) ) - XYZ(1:NOPT) )/(NIMAGE+1)
240:              ENDIF238:              ENDIF
241:    239:    
242: ! bs360: Interpolation using dihedrals.240: ! bs360: Interpolation using dihedrals.
243:              IF (CHRMMT.AND.CHICDNEB.AND.ICINTERPT) THEN241:              IF (CHRMMT.AND.CHICDNEB.AND.ICINTERPT) THEN
244:                 QS(1:NOPT)=XYZ(1:NOPT)         242:                 QS(1:NOPT)=XYZ(1:NOPT)         
245:                 QF(1:NOPT)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2))243:                 QF(1:NOPT)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2))
246:              ENDIF244:              ENDIF
247: ! end bs360245: ! end bs360
248: ! msb50: Interpolation using dihedrals for amber246: 
249:              IF (AMBERT.AND.AMBICDNEBT.AND.AMBERICT) THEN 
250:                 QS(1:NOPT)=XYZ(1:NOPT) 
251:                 QF(1:NOPT)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2)) 
252:                 PRINT*, "msb50 NOPT", NOPT 
253: !                DO K=1, 3*NATOMS 
254: !                    PRINT '(6f9.3)', QF(6*(K-1)+1:6*K) 
255: !                ENDDO 
256:              ENDIF 
257: ! end msb50 
258:              DO I=2,NIMAGE+1247:              DO I=2,NIMAGE+1
259:                   XYZ(NOPT*(I-1)+1:NOPT*I) = XYZ(1:NOPT) + DELTAX*(I-1)248:                   XYZ(NOPT*(I-1)+1:NOPT*I) = XYZ(1:NOPT) + DELTAX*(I-1)
260:        249:        
261: ! bs360250: ! bs360
262:                   IF (CHRMMT.AND.CHICDNEB.AND.ICINTERPT) THEN251:                   IF (CHRMMT.AND.CHICDNEB.AND.ICINTERPT) THEN
263:                      COORDS(1:NOPT)=XYZ(NOPT*(I-1)+1:NOPT*I)252:                      COORDS(1:NOPT)=XYZ(NOPT*(I-1)+1:NOPT*I)
264:                      SFRAC=1.D0-(I-1.D0)/(NIMAGE+1.D0)                  253:                      SFRAC=1.D0-(I-1.D0)/(NIMAGE+1.D0)                  
265: !                     write(6,*) 'i,NIMAGE+1,sfrac= ',i,NIMAGE+1,sfrac254: !                     write(6,*) 'i,NIMAGE+1,sfrac= ',i,NIMAGE+1,sfrac
266:                      CALL ICINTERPOL(COORDS,QS,QF,SFRAC)255:                      CALL ICINTERPOL(COORDS,QS,QF,SFRAC)
267:                      XYZ(NOPT*(I-1)+1:NOPT*I)=COORDS(1:NOPT)        256:                      XYZ(NOPT*(I-1)+1:NOPT*I)=COORDS(1:NOPT)        
268: !                     write(100+I,'(3F10.3)') XYZ(NOPT*(I-1)+1:NOPT*I)257: !                     write(100+I,'(3F10.3)') XYZ(NOPT*(I-1)+1:NOPT*I)
269: !                     call flush(100+I) 
270:                      call flush(100+I,ISTAT)258:                      call flush(100+I,ISTAT)
271:                   ENDIF259:                   ENDIF
272: ! end bs360260: ! end bs360
273:  
274: ! msb50 
275:                   IF (AMBERT.AND.AMBICDNEBT.AND.AMBERICT) THEN 
276:                      COORDS(1:NOPT)=XYZ(NOPT*(I-1)+1:NOPT*I) 
277:                      SFRAC=1.D0-(I-1.D0)/(NIMAGE+1.D0) 
278:                      IF (.NOT.ALLOCATED(NICTOT)) THEN 
279:                          CALL SETDIHEAM() 
280:                      ENDIF 
281:                      CALL TAKESTEPAMDIHED(COORDS, QS, QF,SFRAC) 
282:                      write(6,*) 'i,NIMAGE+1,sfrac= ',i,NIMAGE+1,sfrac 
283:                      XYZ(NOPT*(I-1)+1:NOPT*I)=COORDS(1:NOPT) 
284:                      !DO k=1,(NOPT/3) 
285:                      !  PRINT*, XYZ(3*(k-1)+1:3*k)  
286:                      !ENDDO 
287: !                    write(100+I,'(3F10.3)') XYZ(NOPT*(I-1)+1:NOPT*I) 
288: !                     call flush(100+I) 
289:                   ENDIF 
290: ! end msb50 
291:  
292:                   IF (DESMINT) THEN261:                   IF (DESMINT) THEN
293:                      PREVDIH => DIHINFO(I,:)262:                      PREVDIH => DIHINFO(I,:)
294:                      CALL TRANSBACKDELTA(DELTAX,DELTACART,XYZCART(3*NATOMS*(I-2)+1:3*NATOMS*(I-1)),&263:                      CALL TRANSBACKDELTA(DELTAX,DELTACART,XYZCART(3*NATOMS*(I-2)+1:3*NATOMS*(I-1)),&
295:                           & NINTC,3*NATOMS,NNZ,KD,FAILED,INTPTEST,INTEPSILON)264:                           & NINTC,3*NATOMS,NNZ,KD,FAILED,INTPTEST,INTEPSILON)
296:                      XYZCART(3*NATOMS*(I-1)+1:3*NATOMS*I) = XYZCART(3*NATOMS*(I-2)+1:3*NATOMS*(I-1)) + DELTACART(1:3*NATOMS)265:                      XYZCART(3*NATOMS*(I-1)+1:3*NATOMS*I) = XYZCART(3*NATOMS*(I-2)+1:3*NATOMS*(I-1)) + DELTACART(1:3*NATOMS)
297:                      CALL NEWMINDIST(XYZCART(1:3*NATOMS),XYZCART(3*NATOMS*(I-1)+1:3*NATOMS*I),&266:                      CALL NEWMINDIST(XYZCART(1:3*NATOMS),XYZCART(3*NATOMS*(I-1)+1:3*NATOMS*I),&
298:                           & NATOMS,DIST,.FALSE.,.FALSE.,ZSYM(1),.FALSE.,.FALSE.,.FALSE.,RMAT)267:                           & NATOMS,DIST,.FALSE.,.FALSE.,ZSYM(1),.FALSE.,.FALSE.,.FALSE.,RMAT)
299:                   ENDIF268:                   ENDIF
300:              ENDDO269:              ENDDO
301:              270:              


r15547/output.f90 2017-01-21 10:33:10.882250000 +0000 r15546/output.f90 2017-01-21 10:33:17.674250000 +0000
119: !              dummy=>dummy%next119: !              dummy=>dummy%next
120: !         enddo120: !         enddo
121: !         write(*,'(a)') '.'121: !         write(*,'(a)') '.'
122: 122: 
123:           ! ------ bs360 : more general printout (without MaxPrintOut) -------123:           ! ------ bs360 : more general printout (without MaxPrintOut) -------
124:           WRITE(*,'(1x,a,i4,a)',advance='No') 'Following ',nt,' images are candidates for TS:'124:           WRITE(*,'(1x,a,i4,a)',advance='No') 'Following ',nt,' images are candidates for TS:'
125:           DO J=1,NTSMAX125:           DO J=1,NTSMAX
126:                WRITE(*,'(i5)',advance='No') dummy%i-1126:                WRITE(*,'(i5)',advance='No') dummy%i-1
127:                IF (.NOT.ASSOCIATED(DUMMY%NEXT)) EXIT127:                IF (.NOT.ASSOCIATED(DUMMY%NEXT)) EXIT
128:                DUMMY=>DUMMY%NEXT128:                DUMMY=>DUMMY%NEXT
129:                !msb50 
130:           ENDDO129:           ENDDO
131:  
132:           PRINT *,' '130:           PRINT *,' '
133:           ! ------ end bs360 ---------------------------131:           ! ------ end bs360 ---------------------------
134:           132:           
135:           IF (OPTIMIZETS) THEN133:           IF (OPTIMIZETS) THEN
136:                WRITE(*,'(1x,a)',advance='No') 'Converged to TS (number of iterations):     '134:                WRITE(*,'(1x,a)',advance='No') 'Converged to TS (number of iterations):     '
137: 135: 
138:                DUMMY=>FIRST136:                DUMMY=>FIRST
139:                NTSFOUND=0137:                NTSFOUND=0
140:                CALL MYCPU_TIME(STARTTIME,.FALSE.)138:                CALL MYCPU_TIME(STARTTIME,.FALSE.)
141:                DO J=1,NTSMAX139:                DO J=1,NTSMAX


r15547/potential.f 2017-01-21 10:33:14.006250000 +0000 r15546/potential.f 2017-01-21 10:33:20.810250000 +0000
188:       IF (CONTAINER) CALL RAD(COORDS)188:       IF (CONTAINER) CALL RAD(COORDS)
189: C     IF (CONTAINER) CALL RAD(COORDS,ENERGY,VNEW,GTEST)189: C     IF (CONTAINER) CALL RAD(COORDS,ENERGY,VNEW,GTEST)
190: 190: 
191:       IF (VARIABLES) THEN191:       IF (VARIABLES) THEN
192: C        CALL CTEST(NATOMS, COORDS, VNEW, ENERGY, GTEST, STEST)192: C        CALL CTEST(NATOMS, COORDS, VNEW, ENERGY, GTEST, STEST)
193: C        CALL TWODFUNC(COORDS,VNEW,ENERGY,GTEST,STEST)193: C        CALL TWODFUNC(COORDS,VNEW,ENERGY,GTEST,STEST)
194: C        CALL MB(COORDS,VNEW,ENERGY,GTEST,STEST)194: C        CALL MB(COORDS,VNEW,ENERGY,GTEST,STEST)
195: C         CALL P4DIFF(NATOMS,COORDS,VNEW,ENERGY,PARAM1,GTEST,STEST)195: C         CALL P4DIFF(NATOMS,COORDS,VNEW,ENERGY,PARAM1,GTEST,STEST)
196:          CALL P4DIFF(NATOMS,COORDS,VNEW,ENERGY,PARAM1,PARAM2,GTEST,STEST)196:          CALL P4DIFF(NATOMS,COORDS,VNEW,ENERGY,PARAM1,PARAM2,GTEST,STEST)
197: 197: 
198: C         DIFF=1.0D-3198: !         DIFF=1.0D-3
199: C         PRINT*,'analytic and numerical gradients: NATOMS=',NATOMS199: !         PRINT*,'analytic and numerical gradients: NATOMS=',NATOMS
200: C         DO J1=1,NATOMS200: !         DO J1=1,NATOMS
201: C            COORDS(J1)=COORDS(J1)+DIFF201: !            COORDS(J1)=COORDS(J1)+DIFF
202: C            CALL P4DIFF(NATOMS,COORDS,VPLUS,EPLUS,PARAM1,.FALSE.,.FALSE.)202: !            CALL P4DIFF(NATOMS,COORDS,VPLUS,EPLUS,PARAM1,.FALSE.,.FALSE.)
203: C            COORDS(J1)=COORDS(J1)-2.0D0*DIFF203: !            COORDS(J1)=COORDS(J1)-2.0D0*DIFF
204: C            CALL P4DIFF(NATOMS,COORDS,VMINUS,EMINUS,PARAM1,.FALSE.,.FALSE.)204: !            CALL P4DIFF(NATOMS,COORDS,VMINUS,EMINUS,PARAM1,.FALSE.,.FALSE.)
205: C            COORDS(J1)=COORDS(J1)+DIFF205: !            COORDS(J1)=COORDS(J1)+DIFF
206: C            IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*ABS((VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.0.0D0)) THEN206: !            IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*ABS((VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.0.0D0)) THEN
207: C               WRITE(*,'(I5,2G20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)207: !               WRITE(*,'(I5,2G20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)
208: C            ENDIF208: !            ENDIF
209: C         ENDDO209: !         ENDDO
210: C         PRINT*,'analytic and numerical second derivatives:'210: !         PRINT*,'analytic and numerical second derivatives:'
211: C         DO J1=1,NATOMS211: !         DO J1=1,NATOMS
212: C            COORDS(J1)=COORDS(J1)+DIFF212: !            COORDS(J1)=COORDS(J1)+DIFF
213: C            CALL P4DIFF(NATOMS,COORDS,VPLUS,EPLUS,PARAM1,.TRUE.,.FALSE.)213: !            CALL P4DIFF(NATOMS,COORDS,VPLUS,EPLUS,PARAM1,.TRUE.,.FALSE.)
214: C            COORDS(J1)=COORDS(J1)-2.0D0*DIFF214: !            COORDS(J1)=COORDS(J1)-2.0D0*DIFF
215: C            CALL P4DIFF(NATOMS,COORDS,VMINUS,EMINUS,PARAM1,.TRUE.,.FALSE.)215: !            CALL P4DIFF(NATOMS,COORDS,VMINUS,EMINUS,PARAM1,.TRUE.,.FALSE.)
216: C            COORDS(J1)=COORDS(J1)+DIFF216: !            COORDS(J1)=COORDS(J1)+DIFF
217: C            DO J2=1,NATOMS217: !            DO J2=1,NATOMS
218: C               IF (ABS(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)).GT.1.0D-1) THEN218: !               IF (ABS(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)).GT.1.0D-1) THEN
219: C                  WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'   X'219: !                  WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'   X'
220: C               ELSE220: !               ELSE
221: C                  WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)221: !                  WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)
222: C               ENDIF222: !               ENDIF
223: C            ENDDO223: !            ENDDO
224: C         ENDDO224: !         ENDDO
225: 225: 
226:          IF (PTEST) THEN226:          IF (PTEST) THEN
227:             WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '227:             WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
228:             WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '228:             WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
229:          ENDIF229:          ENDIF
230: 230: 
231: C         IF (RESTART) THEN231: C         IF (RESTART) THEN
232: CC232: CC
233: CC  This is an inefficient dirty fix233: CC  This is an inefficient dirty fix
234: CC234: CC
323:       ELSE IF (SIO2T) THEN323:       ELSE IF (SIO2T) THEN
324:          CALL SIO2(NATOMS, COORDS, VNEW, ENERGY, PARAM1, PARAM2, PARAM3, PARAM4, ZSYM, GTEST, SSTEST)324:          CALL SIO2(NATOMS, COORDS, VNEW, ENERGY, PARAM1, PARAM2, PARAM3, PARAM4, ZSYM, GTEST, SSTEST)
325:          IF (SIO2C6T) THEN325:          IF (SIO2C6T) THEN
326:             CALL SIO2C6(NATOMS, COORDS, VNEW, EDISP, C6PP, C6MM, C6PM, ZSYM, GTEST, SSTEST)326:             CALL SIO2C6(NATOMS, COORDS, VNEW, EDISP, C6PP, C6MM, C6PM, ZSYM, GTEST, SSTEST)
327:             IF (PTEST) WRITE(*,'(A,F20.10,A)') ' Dispersion energy=',EDISP,' hartree'327:             IF (PTEST) WRITE(*,'(A,F20.10,A)') ' Dispersion energy=',EDISP,' hartree'
328:             ENERGY=ENERGY+EDISP328:             ENERGY=ENERGY+EDISP
329:          ENDIF329:          ENDIF
330: C330: C
331: C  Check numerical first and second derivatives331: C  Check numerical first and second derivatives
332: C332: C
333: C        DIFF=1.0D-4333: !        DIFF=1.0D-4
334: C        PRINT*,'analytic and numerical gradients:'334: !        PRINT*,'analytic and numerical gradients:'
335: C        DO J1=1,3*NATOMS335: !        DO J1=1,3*NATOMS
336: C           COORDS(J1)=COORDS(J1)+DIFF336: !           COORDS(J1)=COORDS(J1)+DIFF
337: C           CALL SIO2(NATOMS,COORDS,VPLUS,EPLUS,PARAM1, PARAM2, PARAM3, PARAM4,ZSYM,.FALSE.,.FALSE.)337: !           CALL SIO2(NATOMS,COORDS,VPLUS,EPLUS,PARAM1, PARAM2, PARAM3, PARAM4,ZSYM,.FALSE.,.FALSE.)
338: C           IF (SIO2C6T) THEN338: !           IF (SIO2C6T) THEN
339: C              CALL SIO2C6(NATOMS, COORDS, VPLUS, EDISP, C6PP, C6MM, C6PM, ZSYM, .FALSE.,.FALSE.)339: !              CALL SIO2C6(NATOMS, COORDS, VPLUS, EDISP, C6PP, C6MM, C6PM, ZSYM, .FALSE.,.FALSE.)
340: C              EPLUS=EPLUS+EDISP340: !              EPLUS=EPLUS+EDISP
341: C           ENDIF341: !           ENDIF
342: C           COORDS(J1)=COORDS(J1)-2.0D0*DIFF342: !           COORDS(J1)=COORDS(J1)-2.0D0*DIFF
343: C           CALL SIO2(NATOMS,COORDS,VMINUS,EMINUS,PARAM1, PARAM2, PARAM3, PARAM4,ZSYM,.FALSE.,.FALSE.)343: !           CALL SIO2(NATOMS,COORDS,VMINUS,EMINUS,PARAM1, PARAM2, PARAM3, PARAM4,ZSYM,.FALSE.,.FALSE.)
344: C           IF (SIO2C6T) THEN344: !           IF (SIO2C6T) THEN
345: C              CALL SIO2C6(NATOMS, COORDS, VMINUS, EDISP, C6PP, C6MM, C6PM, ZSYM, .FALSE.,.FALSE.)345: !              CALL SIO2C6(NATOMS, COORDS, VMINUS, EDISP, C6PP, C6MM, C6PM, ZSYM, .FALSE.,.FALSE.)
346: C              EMINUS=EMINUS+EDISP346: !              EMINUS=EMINUS+EDISP
347: C           ENDIF347: !           ENDIF
348: C           COORDS(J1)=COORDS(J1)+DIFF348: !           COORDS(J1)=COORDS(J1)+DIFF
349: C           IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*(VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1).GT.1.0D0)) THEN349: !           IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*(VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1).GT.1.0D0)) THEN
350: C              WRITE(*,'(I5,2F20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)350: !              WRITE(*,'(I5,2F20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)
351: C           ENDIF351: !           ENDIF
352: C        ENDDO352: !        ENDDO
353: !        PRINT*,'analytic and numerical second derivatives:'353: !        PRINT*,'analytic and numerical second derivatives:'
354: !        DO J1=1,3*NATOMS354: !        DO J1=1,3*NATOMS
355: !           COORDS(J1)=COORDS(J1)+DIFF355: !           COORDS(J1)=COORDS(J1)+DIFF
356: !           CALL SIO2(NATOMS,COORDS,VPLUS,EPLUS,PARAM1, PARAM2, PARAM3, PARAM4,ZSYM,.TRUE.,.FALSE.)356: !           CALL SIO2(NATOMS,COORDS,VPLUS,EPLUS,PARAM1, PARAM2, PARAM3, PARAM4,ZSYM,.TRUE.,.FALSE.)
357: !           IF (SIO2C6T) CALL SIO2C6(NATOMS, COORDS, VPLUS, EDISP, C6PP, C6MM, C6PM, ZSYM, .TRUE.,.FALSE.)357: !           IF (SIO2C6T) CALL SIO2C6(NATOMS, COORDS, VPLUS, EDISP, C6PP, C6MM, C6PM, ZSYM, .TRUE.,.FALSE.)
358: !           COORDS(J1)=COORDS(J1)-2.0D0*DIFF358: !           COORDS(J1)=COORDS(J1)-2.0D0*DIFF
359: !           CALL SIO2(NATOMS,COORDS,VMINUS,EMINUS,PARAM1, PARAM2, PARAM3, PARAM4,ZSYM,.TRUE.,.FALSE.)359: !           CALL SIO2(NATOMS,COORDS,VMINUS,EMINUS,PARAM1, PARAM2, PARAM3, PARAM4,ZSYM,.TRUE.,.FALSE.)
360: !           IF (SIO2C6T) CALL SIO2C6(NATOMS, COORDS, VMINUS, EDISP, C6PP, C6MM, C6PM, ZSYM, .TRUE.,.FALSE.)360: !           IF (SIO2C6T) CALL SIO2C6(NATOMS, COORDS, VMINUS, EDISP, C6PP, C6MM, C6PM, ZSYM, .TRUE.,.FALSE.)
361: !           COORDS(J1)=COORDS(J1)+DIFF361: !           COORDS(J1)=COORDS(J1)+DIFF
362: !           DO J2=1,3*NATOMS362: !           DO J2=1,3*NATOMS
702:          CALL MPDIFF(NATOMS,COORDS,VNEW,ENERGY,PARAM1,PARAM2,PARAM3,PARAM4,PARAM5,GTEST,SSTEST)702:          CALL MPDIFF(NATOMS,COORDS,VNEW,ENERGY,PARAM1,PARAM2,PARAM3,PARAM4,PARAM5,GTEST,SSTEST)
703:          IF (PTEST) WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' epsilon'703:          IF (PTEST) WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' epsilon'
704:          WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' epsilon'704:          WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' epsilon'
705:       ELSE IF (ZSYM(NATOMS).EQ.'GP') THEN705:       ELSE IF (ZSYM(NATOMS).EQ.'GP') THEN
706:          CALL  GUPTA(NATOMS,COORDS,VNEW,ENERGY,GTEST,GUPTATYPE)706:          CALL  GUPTA(NATOMS,COORDS,VNEW,ENERGY,GTEST,GUPTATYPE)
707:          IF (PTEST) WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' eV'707:          IF (PTEST) WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' eV'
708:          WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' eV'708:          WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' eV'
709:       ELSE IF (ZSYM(NATOMS).EQ.'DS') THEN709:       ELSE IF (ZSYM(NATOMS).EQ.'DS') THEN
710:          CALL MPDIFFDS(NATOMS,COORDS,VNEW,ENERGY,PARAM1,PARAM2,PARAM3,PARAM4,PARAM5,GTEST,SSTEST)710:          CALL MPDIFFDS(NATOMS,COORDS,VNEW,ENERGY,PARAM1,PARAM2,PARAM3,PARAM4,PARAM5,GTEST,SSTEST)
711: 711: 
712: C        DIFF=1.0D-5712: !        DIFF=1.0D-5
713: C        PRINT*,'analytic and numerical gradients:'713: !        PRINT*,'analytic and numerical gradients:'
714: C        DO J1=1,3*NATOMS714: !        DO J1=1,3*NATOMS
715: C        DO J1=18,18715: !        DO J1=18,18
716: C           IF (FROZEN((J1-1)/3+1)) CYCLE716: !           IF (FROZEN((J1-1)/3+1)) CYCLE
717: C           COORDS(J1)=COORDS(J1)+DIFF717: !           COORDS(J1)=COORDS(J1)+DIFF
718: C           CALL MPDIFFDS(NATOMS,COORDS,VPLUS,EPLUS,PARAM1,PARAM2,PARAM3,PARAM4,PARAM5,.FALSE.,.FALSE.)718: !           CALL MPDIFFDS(NATOMS,COORDS,VPLUS,EPLUS,PARAM1,PARAM2,PARAM3,PARAM4,PARAM5,.FALSE.,.FALSE.)
719: C           COORDS(J1)=COORDS(J1)-2.0D0*DIFF719: !           COORDS(J1)=COORDS(J1)-2.0D0*DIFF
720: C           CALL MPDIFFDS(NATOMS,COORDS,VMINUS,EMINUS,PARAM1,PARAM2,PARAM3,PARAM4,PARAM5,.FALSE.,.FALSE.)720: !           CALL MPDIFFDS(NATOMS,COORDS,VMINUS,EMINUS,PARAM1,PARAM2,PARAM3,PARAM4,PARAM5,.FALSE.,.FALSE.)
721: C           COORDS(J1)=COORDS(J1)+DIFF721: !           COORDS(J1)=COORDS(J1)+DIFF
722: C           IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*(VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1).GT.1.0D0)) THEN722: !           IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*(VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1).GT.1.0D0)) THEN
723: C              WRITE(*,'(I5,2G20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)723: !              WRITE(*,'(I5,2G20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)
724: C           ENDIF724: !           ENDIF
725: C        ENDDO725: !        ENDDO
726: !        STOP726: !        STOP
727: 727: 
728:          IF (PTEST) WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' epsilon'728:          IF (PTEST) WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' epsilon'
729:          WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' epsilon'729:          WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' epsilon'
730:       ELSE IF (ZSYM(NATOMS).EQ.'MS') THEN730:       ELSE IF (ZSYM(NATOMS).EQ.'MS') THEN
731:          IF (PRESSURE) THEN731:          IF (PRESSURE) THEN
732:             CALL MSLATMIN(NATOMS,COORDS,PARAM1,PARAM2,PARAM3,PARAM4)732:             CALL MSLATMIN(NATOMS,COORDS,PARAM1,PARAM2,PARAM3,PARAM4)
733:             PRINT*,'Lattice constant optimised'733:             PRINT*,'Lattice constant optimised'
734:             PRINT*,'New Box length in x=',PARAM2734:             PRINT*,'New Box length in x=',PARAM2
735:             PRINT*,'New Box length in y=',PARAM3735:             PRINT*,'New Box length in y=',PARAM3
980:          IF (PTEST) THEN980:          IF (PTEST) THEN
981:             WRITE(*,10) ' potential> Energy for last cycle=',ENERGY981:             WRITE(*,10) ' potential> Energy for last cycle=',ENERGY
982:             WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY982:             WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY
983:          ENDIF983:          ENDIF
984:       ELSE IF (ZSYM(NATOMS).EQ.'CD') THEN984:       ELSE IF (ZSYM(NATOMS).EQ.'CD') THEN
985:          CALL FCAPSID(NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST)985:          CALL FCAPSID(NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST)
986:          IF (PTEST) THEN986:          IF (PTEST) THEN
987:             WRITE(*,10) ' potential> Energy for last cycle=',ENERGY987:             WRITE(*,10) ' potential> Energy for last cycle=',ENERGY
988:             WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY988:             WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY
989:          ENDIF989:          ENDIF
990:  
991: C990: C
992: C  Check numerical first and second derivatives991: C  Check numerical first and second derivatives
993: C992: C
994: C        DIFF=1.0D-4993: C        DIFF=1.0D-4
995: C        PRINT*,'analytic and numerical gradients:'994: C        PRINT*,'analytic and numerical gradients:'
996: C        DO J1=1,3*NATOMS995: C        DO J1=1,3*NATOMS
997: C           COORDS(J1)=COORDS(J1)+DIFF996: C           COORDS(J1)=COORDS(J1)+DIFF
998: C           CALL FCAPSID(NATOMS,COORDS,HDUMM,VPLUS,EPLUS,.FALSE.,.FALSE.)997: C           CALL FCAPSID(NATOMS,COORDS,HDUMM,VPLUS,EPLUS,.FALSE.,.FALSE.)
999: C           COORDS(J1)=COORDS(J1)-2.0D0*DIFF998: C           COORDS(J1)=COORDS(J1)-2.0D0*DIFF
1000: C           CALL FCAPSID(NATOMS,COORDS,HDUMM,VPLUS,EMINUS,.FALSE.,.FALSE.)999: C           CALL FCAPSID(NATOMS,COORDS,HDUMM,VPLUS,EMINUS,.FALSE.,.FALSE.)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0