hdiff output

r28580/amberinterface.f 2017-03-30 13:30:32.312543027 +0100 r28579/amberinterface.f 2017-03-30 13:30:34.100566797 +0100
1414: if (.not.fusion) exit1414: if (.not.fusion) exit
1415: ! write(222,'(3F20.10)') x1,y1,z1 1415: ! write(222,'(3F20.10)') x1,y1,z1 
1416: END DO1416: END DO
1417: ! write(MYUNITNEW,*) 'repositioned some atoms, nsteps=', k1417: ! write(MYUNITNEW,*) 'repositioned some atoms, nsteps=', k
1418: !close(222)1418: !close(222)
1419: !STOP1419: !STOP
1420: end subroutine checkdistances1420: end subroutine checkdistances
1421: 1421: 
1422: subroutine check_cistrans_rna(coords,natoms,atomlabels,goodstructure)1422: subroutine check_cistrans_rna(coords,natoms,atomlabels,goodstructure)
1423: use modamber9,only : myunitnew1423: use modamber9,only : myunitnew
1424:  
1425: implicit none1424: implicit none
1426: 1425: 
1427: integer,intent(in) :: natoms1426: integer,intent(in) :: natoms
1428: character(len=5),intent(in) :: atomlabels(natoms)1427: character(len=5),intent(in) :: atomlabels(natoms)
1429: integer :: i,j,l1428: integer :: i,j
1430: double precision,intent(in) :: coords(3*natoms)1429: double precision,intent(in) :: coords(3*natoms)
1431: double precision, dimension(1,3) :: vector1, vector2 
1432: double precision, dimension(3) :: cross_prod 
1433: double precision :: moduli 
1434: double precision, parameter :: pi = 4.0*atan(1.0) 
1435: double precision, dimension(15,3) :: b_vector 
1436: double precision, dimension(15,3) :: v4 
1437: ! stuff for detecting cis-trans isomerisation in ribose rings 
1438: integer :: resid  
1439: integer, dimension(20) :: atom 
1440: integer, dimension(5,4) :: dihed 
1441: double precision, dimension(15) :: distance 
1442: double precision :: angle(3) 
1443: logical, intent(out) :: goodstructure 
1444: double precision, dimension(5) :: dihedral_angle 
1445: 1430: 
 1431: 
 1432: 
 1433: ! stuff for detecting cis-trans isomerisation in ribose rings
 1434: integer :: dihed(3,4),resid, atom(9)
 1435: double precision :: angle(3),distance(14)
 1436: logical,intent(out) :: goodstructure
1446: goodstructure=.TRUE.1437: goodstructure=.TRUE.
1447: 1438: 
1448:       resid=01439:       resid=0
1449:       do i=1,natoms   1440:       do i=1,natoms   
1450:  
1451: !THE ALTERNATE CHIRAL CHECKS ARE TURNED OFF FOR NOW           
1452: ! new cis-trans checks based on torsion angles of the chiral carbons. 
1453: ! The torsions to be checked are according to the definition from Yildirim 
1454: ! et.al., J.Phys.Chem B, 2010. 
1455: ! C2'-O4'-N9/N1-H1' 
1456: ! C3'-C1'-H2'1-O2' 
1457: ! O3'-C2'-C4'-H3' 
1458: ! C3'-C5'-O4'-H4' 
1459: ! C4'-O5'-H5'1-H5'2 
1460: ! The arctangents of the  dihedral angles should always be positive 
1461: ! The range is chosen initially from the distribution sampled in TIP3P/MD.  
1462: ! If there is a sign change then it means that there is an undesired 
1463: ! chiral inversion. 
1464: ! select atoms the distances of which we are going to check 1441: ! select atoms the distances of which we are going to check 
1465:  
1466: !        if (atomlabels(i).eq."O5'")then 
1467: !         atom(18) = i   
1468: !         elseif (atomlabels(i).eq."C5'")then 
1469: !         atom(14) = i 
1470: !         elseif (atomlabels(i).eq. "H5'1")then 
1471: !         atom(19) = i 
1472: !         elseif (atomlabels(i).eq. "H5'2")then 
1473: !         atom(20) = i 
1474: !         elseif (atomlabels(i).eq. "C4'")then 
1475: !         atom(11) = i 
1476:          !counter = counter + 1 
1477: !         atom(17) = i 
1478: !         elseif (atomlabels(i).eq."H4'")then 
1479: !         atom(16) = i 
1480: !         elseif (atomlabels(i).eq. "O4'")then 
1481: !         atom(2) = i 
1482: !         atom(15) = i 
1483: !         elseif (atomlabels(i).eq. "C1'")then 
1484: !         atom(6) = i 
1485: !         elseif (atomlabels(i).eq."H1'")then 
1486: !         atom(4) = i 
1487: !         atom(3) = i+1 
1488: !         elseif (atomlabels(i).eq."C3'")then 
1489: !         atom(5) = i 
1490:          !counter = counter + 1 
1491: !         atom(13) = i 
1492:          !counter = counter + 1 
1493:          !dihed(4,2) = i 
1494:          !endif     
1495: !         elseif (atomlabels(i).eq. "H3'")then 
1496: !         atom(12) = i 
1497: !         elseif (atomlabels(i).eq."C2'")then 
1498: !         atom(1) = i  
1499: !         atom(10) = i 
1500: !         elseif (atomlabels(i).eq. "H2'1")then 
1501: !         atom(7) = i 
1502: !         elseif (atomlabels(i).eq. "O2'")then 
1503: !         atom(8) = i 
1504: !         elseif (atomlabels(i).eq. "O3'")then 
1505: !         atom(9) = i 
1506: !         resid = resid + 1 
1507: !Now compute the bond vectors and the relevant unit vectors 
1508: !And compute the dihedrals using the atan2 function as it returns angle 
1509: !With the correct sign 
1510: ! Calculates the required unit bond vectors required for the subsequent 
1511: ! cross products  
1512: ! cross products between b1 and b2; b2 and b3; b2 and v6 
1513: ! where v6 is the normalized cross product between b1 and b2 
1514: ! the cross-products also need to be normalized. 
1515: ! the cross-product subroutine should return cross-prod with the normalizations. 
1516: !   DO l = 1,5 
1517: !      call bond_vector(atom(l*4-3),atom(l*4-2),coords,natoms,b_vector(l*3-2,1),b_vector(l*3-2,2),b_vector(l*3-2,3)) 
1518: !      call bond_vector(atom(l*4-2),atom(l*4-1),coords,natoms,b_vector(l*3-1,1),b_vector(l*3-1,2),b_vector(l*3-1,3)) 
1519: !      call bond_vector(atom(l*4-1),atom(l*4),coords,natoms,b_vector(l*3,1),b_vector(l*3,2),b_vector(l*3,3)) 
1520: !   enddo 
1521:           
1522: ! DO THE CROSS PRODUCTS 
1523: !DO l = 1,5 
1524: ! call cross_prod_na(b_vector(l*3-2,:),b_vector(l*3-1,:),cross_prod,moduli)  
1525: ! v4(l*3-2,1) = cross_prod(1)/moduli 
1526: ! v4(l*3-2,2) = cross_prod(2)/moduli 
1527: ! v4(l*3-2,3) = cross_prod(3)/moduli 
1528: !call cross_prod_na(b_vector(l*3-1,:),b_vector(l*3,:),cross_prod,moduli) 
1529: ! v4(l*3-1,1) = cross_prod(1)/moduli 
1530: ! v4(l*3-1,2) = cross_prod(2)/moduli 
1531: !  v4(l*3-1,3) = cross_prod(3)/moduli 
1532: !call cross_prod_na(b_vector(l*3-1,:),v4(l*3-2,:),cross_prod,moduli)  
1533: ! v4(l*3,1) = cross_prod(1)/moduli 
1534: ! v4(l*3,2) = cross_prod(2)/moduli 
1535: ! v4(l*3,3) = cross_prod(3)/moduli  
1536: !dihedral_angle(l) = atan2(DOT_PRODUCT(v4(l*3,:),v4(l*3-1,:)),DOT_PRODUCT(v4(l*3-2,:),v4(l*3-1,:))) * (180/pi) 
1537: !ENDDO    
1538: ! ALLOWING FOR METHYLENE GROUPS IN THE BACKBONE TO ROTATE. IF THEY HAVE TO BE 
1539: ! STOPPED PUT l = 5.  
1540:  
1541: ! DO l = 1,5 
1542: ! IF(dihedral_angle(l).le.0)THEN 
1543: ! WRITE(MYUNITNEW,'(A,I3,A,I3,A,F8.3)')"chiralcheck_rna> Sign change for angle residue ",resid," angle ",l," value ", dihedral_angle(l) 
1544: ! goodstructure=.false. 
1545: ! ENDIF 
1546: ! ENDDO 
1547:  
1548:  
1549:       
1550:  !if(dihedral_angle(j).lt.0)then 
1551:   !  goodstructure=.false. 
1552:    ! endif 
1553: !WRITE(MYUNITNEW,'(A,I3,A,I3,A,F8.3)') "newchiralcheck_rna> dihedral angle  residue",resid,"is angle",j,"value",dihedral_angle(j) 
1554:  !        endif 
1555:  !      enddo 
1556: !if(goodstructure)then 
1557: !write(MYUNITNEW,'(A)')"structure passed the torsion tests"   
1558: !endif 
1559:    !Commenting out the lines below just for checking ? 
1560: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
1561: ! REVERTING BACK TO SZILARD'S ORIGINAL CHECKS                   
1562:  
1563:         if      (atomlabels(i).eq."C5'") then1442:         if      (atomlabels(i).eq."C5'") then
1564:                 atom(1)=i1443:                 atom(1)=i
1565:         else if (atomlabels(i).eq."H4'") then 1444:         else if (atomlabels(i).eq."H4'") then 
1566:                 dihed(1,1)=i1445:                 dihed(1,1)=i
1567:                 atom(2)=i1446:                 atom(2)=i
1568:         else if (atomlabels(i).eq."N9".or.atomlabels(i).eq."N1") then 1447:         else if (atomlabels(i).eq."N9".or.atomlabels(i).eq."N1") then 
1569:                 atom(9)=i1448:                 atom(9)=i
1570:         else if (atomlabels(i).eq."C4'") then1449:         else if (atomlabels(i).eq."C4'") then
1571:                 dihed(1,2)=i1450:                 dihed(1,2)=i
1572:         else if (atomlabels(i).eq."C3'") then1451:         else if (atomlabels(i).eq."C3'") then
1576:                 dihed(1,4)=i                    ! H4'-C4'-C3'-H3' should be trans1455:                 dihed(1,4)=i                    ! H4'-C4'-C3'-H3' should be trans
1577:                         dihed(2,1)=i1456:                         dihed(2,1)=i
1578:                 atom(3)=i1457:                 atom(3)=i
1579:         else if (atomlabels(i).eq."C2'") then1458:         else if (atomlabels(i).eq."C2'") then
1580:                         dihed(2,3)=i1459:                         dihed(2,3)=i
1581:                                 dihed(3,2)=i1460:                                 dihed(3,2)=i
1582:         else if (atomlabels(i).eq."H2'1") then1461:         else if (atomlabels(i).eq."H2'1") then
1583:                         dihed(2,4)=i            ! H3'-C3'-C2'-H2' should be cis1462:                         dihed(2,4)=i            ! H3'-C3'-C2'-H2' should be cis
1584:                                 dihed(3,1)=i1463:                                 dihed(3,1)=i
1585:                 atom(5)=i1464:                 atom(5)=i
1586:        else if (atomlabels(i).eq."C1'") then1465:         else if (atomlabels(i).eq."C1'") then
1587:                                 dihed(3,3)=i1466:                                 dihed(3,3)=i
1588:         else if (atomlabels(i).eq."H1'") then1467:         else if (atomlabels(i).eq."H1'") then
1589:                                 dihed(3,4)=i    ! H2'-C2'-C1'-H1' should be trans1468:                                 dihed(3,4)=i    ! H2'-C2'-C1'-H1' should be trans
1590:                 atom(8)=i1469:                 atom(8)=i
1591:                 atom(7)=i+1 ! the N of the base follows the H1' atom in the topology file1470:                 atom(7)=i+1 ! the N of the base follows the H1' atom in the topology file
1592:         else if (atomlabels(i).eq."O2'") then1471:         else if (atomlabels(i).eq."O2'") then
1593:                 atom(6)=i1472:                 atom(6)=i
1594:         else if (atomlabels(i).eq."O3'") then1473:         else if (atomlabels(i).eq."O3'") then
1595:                 atom(4)=i1474:                 atom(4)=i
1596:                 resid=resid+11475:                 resid=resid+1
1620: 1499: 
1621:                 end if1500:                 end if
1622:         end do1501:         end do
1623: !        if(abs(angle(1))<80.0D0 .or. abs(angle(2))>90.0D0 .or. abs(angle(3))<80.0D0) then1502: !        if(abs(angle(1))<80.0D0 .or. abs(angle(2))>90.0D0 .or. abs(angle(3))<80.0D0) then
1624: !            write(MYUNITNEW,'(A,I6,A,I3,A,3F8.3)') 'minimum ',k,' residue ',resid,' dihedrals: ',angle(1:3)1503: !            write(MYUNITNEW,'(A,I6,A,I3,A,3F8.3)') 'minimum ',k,' residue ',resid,' dihedrals: ',angle(1:3)
1625: !            goodstructure(k)=.false.1504: !            goodstructure(k)=.false.
1626: !        end if1505: !        end if
1627:        end if1506:        end if
1628:       end do1507:       end do
1629:       if (goodstructure) then1508:       if (goodstructure) then
1630:          write(MYUNITNEW,'(A)') 'structure OK'1509:           write(MYUNITNEW,'(A)') 'structure OK'
1631:      end if1510:       end if
1632: 1511: 
1633: end subroutine check_cistrans_rna1512: end subroutine check_cistrans_rna
1634: 1513: 
1635: SUBROUTINE bond_vector(i,j,coords,natoms,vector_x,vector_y,vector_z) 
1636: IMPLICIT NONE 
1637: !INTEGER, DIMENSION(20),INTENT(IN) :: atom 
1638: !DOUBLE PRECISION, DIMENSION(15,3) :: b_vector 
1639: !DOUBLE PRECISION, DIMENSION(15) :: distance 
1640: INTEGER, INTENT(IN) :: i,j 
1641: !INTEGER :: i 
1642: INTEGER,INTENT(IN) :: natoms 
1643: DOUBLE PRECISION, DIMENSION(3*natoms), INTENT(IN) :: coords 
1644: DOUBLE PRECISION :: vector_x,vector_y,vector_z 
1645: DOUBLE PRECISION :: distance_x, distance_y, distance_z 
1646: distance_x = sqrt((coords(3*j-2)-coords(3*i-2))**2 + (coords(3*j-1)-coords(3*i-1))**2 + (coords(3*j)-coords(3*i))**2) 
1647: vector_x = (coords(3*j-2)-coords(3*i-2))/distance_x 
1648: vector_y = (coords(3*j-1)-coords(3*i-1))/distance_x 
1649: vector_z = (coords(3*j) - coords(3*i))/distance_x 
1650: !WRITE(*,*) distance_x 
1651: END SUBROUTINE bond_vector 
1652:  
1653: subroutine cross_prod_na(vector1,vector2,cross_prod,moduli) 
1654: implicit none 
1655: double precision, dimension(1,3),intent(in) :: vector1,vector2  
1656: ! The first dimension identifies the bond 
1657: double precision, dimension(3),intent(out) :: cross_prod 
1658: double precision, intent(out) :: moduli 
1659:  
1660: cross_prod(1) = vector1(1,2)*vector2(1,3) - vector1(1,3)*vector2(1,2)  
1661: cross_prod(2) = vector1(1,3)*vector2(1,1) - vector1(1,1)*vector2(1,3) 
1662: cross_prod(3) = vector1(1,1)*vector2(1,2) - vector2(1,1)*vector1(1,2) 
1663: moduli = sqrt((cross_prod(1)*cross_prod(1)) + (cross_prod(2)*cross_prod(2)) + & 
1664:          &  (cross_prod(3)*cross_prod(3))) 
1665:  
1666: end subroutine cross_prod_na 
1667:  
1668:   
1669:   
1670:  
1671:  
1672:  
1673:  
1674:  
1675: subroutine check_cistrans_dna(coords,natoms,atomlabels,goodstructure)1514: subroutine check_cistrans_dna(coords,natoms,atomlabels,goodstructure)
1676: use modamber9,only : myunitnew1515: use modamber9,only : myunitnew
1677: implicit none1516: implicit none
 1517: 
1678: integer,intent(in) :: natoms1518: integer,intent(in) :: natoms
1679: character(len=5),intent(in) :: atomlabels(natoms)1519: character(len=5),intent(in) :: atomlabels(natoms)
1680: integer :: i,j,l1520: integer :: i,j
1681: double precision,intent(in) :: coords(3*natoms)1521: double precision,intent(in) :: coords(3*natoms)
1682: double precision, dimension(1,3) :: vector1, vector2 
1683: double precision, dimension(3) :: cross_prod 
1684: double precision :: moduli 
1685: double precision, parameter :: pi = 4.0*atan(1.0) 
1686: double precision, dimension(15,3) :: b_vector 
1687: double precision, dimension(15,3) :: v4 
1688: ! stuff for detecting cis-trans isomerisation in deoxy-ribose rings 
1689: integer :: resid  
1690: integer, dimension(20) :: atom 
1691: integer, dimension(5,4) :: dihed 
1692: double precision, dimension(15) :: distance 
1693: double precision :: angle(3) 
1694: logical, intent(out) :: goodstructure 
1695: double precision, dimension(5) :: dihedral_angle 
1696:  
1697: !integer,intent(in) :: natoms 
1698: !character(len=5),intent(in) :: atomlabels(natoms) 
1699: !integer :: i,j 
1700: !double precision,intent(in) :: coords(3*natoms) 
1701: 1522: 
1702: 1523: 
1703: 1524: 
1704: ! stuff for detecting cis-trans isomerisation in ribose rings1525: ! stuff for detecting cis-trans isomerisation in ribose rings
1705: !integer :: dihed(3,4),resid, atom(10)1526: integer :: dihed(3,4),resid, atom(10)
1706: !double precision :: angle(3),distance(12)1527: double precision :: angle(3),distance(12)
1707: !logical,intent(out) :: goodstructure1528: logical,intent(out) :: goodstructure
1708: goodstructure=.TRUE.1529: goodstructure=.TRUE.
1709: 1530: 
1710:      resid=01531:      resid=0
1711:       do i=1,natoms1532:       do i=1,natoms   
1712: ! if (atomlabels(i).eq."O5'")then 
1713:  !        atom(18) = i   
1714:   !       elseif (atomlabels(i).eq."C5'")then 
1715:    !      atom(14) = i 
1716:     !     elseif (atomlabels(i).eq. "H5'1")then 
1717:    !      atom(19) = i 
1718:    !      elseif (atomlabels(i).eq. "H5'2")then 
1719:    !      atom(20) = i 
1720:    !      elseif (atomlabels(i).eq. "C4'")then 
1721:    !      atom(11) = i 
1722:          !counter = counter + 1 
1723:     !     atom(17) = i 
1724:     !     elseif (atomlabels(i).eq."H4'")then 
1725:     !     atom(16) = i 
1726:     !     elseif (atomlabels(i).eq. "O4'")then 
1727:     !     atom(2) = i 
1728:     !     atom(15) = i 
1729:     !     elseif (atomlabels(i).eq. "C1'")then 
1730:     !     atom(6) = i 
1731:     !     elseif (atomlabels(i).eq."H1'")then 
1732:     !     atom(4) = i 
1733:     !     atom(3) = i+1 
1734:     !     elseif (atomlabels(i).eq."C3'")then 
1735:     !     atom(5) = i 
1736:     !     !counter = counter + 1 
1737:     !     atom(13) = i 
1738:          !counter = counter + 1 
1739:          !dihed(4,2) = i 
1740:          !endif     
1741:     !     elseif (atomlabels(i).eq. "H3'")then 
1742:     !     atom(12) = i 
1743:     !     elseif (atomlabels(i).eq."C2'")then 
1744:     !     atom(1) = i  
1745:     !     atom(10) = i 
1746:     !     elseif (atomlabels(i).eq. "H2'1")then 
1747:     !     atom(7) = i 
1748:     !     elseif (atomlabels(i).eq. "H2'2")then 
1749:     !     atom(8) = i 
1750:     !     elseif (atomlabels(i).eq. "O3'")then 
1751:     !     atom(9) = i 
1752:     !     resid = resid + 1    
1753:  
1754: !DO l = 1,5 
1755:  !     call bond_vector(atom(l*4-3),atom(l*4-2),coords,natoms,b_vector(l*3-2,1),b_vector(l*3-2,2),b_vector(l*3-2,3)) 
1756:   !    call bond_vector(atom(l*4-2),atom(l*4-1),coords,natoms,b_vector(l*3-1,1),b_vector(l*3-1,2),b_vector(l*3-1,3)) 
1757:    !   call bond_vector(atom(l*4-1),atom(l*4),coords,natoms,b_vector(l*3,1),b_vector(l*3,2),b_vector(l*3,3)) 
1758:  !  enddo 
1759:           
1760: ! DO THE CROSS PRODUCTS 
1761: !DO l = 1,5 
1762: ! call cross_prod_na(b_vector(l*3-2,:),b_vector(l*3-1,:),cross_prod,moduli)  
1763: ! v4(l*3-2,1) = cross_prod(1)/moduli 
1764: ! v4(l*3-2,2) = cross_prod(2)/moduli 
1765: ! v4(l*3-2,3) = cross_prod(3)/moduli 
1766: !call cross_prod_na(b_vector(l*3-1,:),b_vector(l*3,:),cross_prod,moduli) 
1767: ! v4(l*3-1,1) = cross_prod(1)/moduli 
1768: ! v4(l*3-1,2) = cross_prod(2)/moduli 
1769: !  v4(l*3-1,3) = cross_prod(3)/moduli 
1770: !call cross_prod_na(b_vector(l*3-1,:),v4(l*3-2,:),cross_prod,moduli)  
1771: ! v4(l*3,1) = cross_prod(1)/moduli 
1772: ! v4(l*3,2) = cross_prod(2)/moduli 
1773: ! v4(l*3,3) = cross_prod(3)/moduli  
1774: !dihedral_angle(l) = atan2(DOT_PRODUCT(v4(l*3,:),v4(l*3-1,:)),DOT_PRODUCT(v4(l*3-2,:),v4(l*3-1,:))) * (180/pi) 
1775: !ENDDO    
1776: ! DO l = 1,4 
1777: ! IF(dihedral_angle(l).le.0)THEN 
1778: ! WRITE(MYUNITNEW,'(A,I3,A,I3,A,F8.3)')"chiralcheck_dna> Sign change in angle for residue ",resid," angle ",l," value ", dihedral_angle(l) 
1779: ! goodstructure=.false. 
1780: ! ENDIF 
1781:  !ENDDO 
1782:  
1783: ! select atoms the distances of which we are going to check 1533: ! select atoms the distances of which we are going to check 
1784:         if      (atomlabels(i).eq."C5'") then1534:         if      (atomlabels(i).eq."C5'") then
1785:                 atom(1)=i1535:                 atom(1)=i
1786:         else if (atomlabels(i).eq."H4'") then 1536:         else if (atomlabels(i).eq."H4'") then 
1787:                 dihed(1,1)=i1537:                 dihed(1,1)=i
1788:                 atom(2)=i1538:                 atom(2)=i
1789:         else if (atomlabels(i).eq."N9") then 1539:         else if (atomlabels(i).eq."N9") then 
1790:                 atom(9)=i1540:                 atom(9)=i
1791:         else if (atomlabels(i).eq."N1") then 1541:         else if (atomlabels(i).eq."N1") then 
1792:                 atom(10)=i1542:                 atom(10)=i
1797:                         dihed(2,2)=i1547:                         dihed(2,2)=i
1798:         else if (atomlabels(i).eq."H3'") then1548:         else if (atomlabels(i).eq."H3'") then
1799:                 dihed(1,4)=i                    ! H4'-C4'-C3'-H3' should be trans1549:                 dihed(1,4)=i                    ! H4'-C4'-C3'-H3' should be trans
1800:                         dihed(2,1)=i1550:                         dihed(2,1)=i
1801:                 atom(3)=i1551:                 atom(3)=i
1802:         else if (atomlabels(i).eq."C2'") then1552:         else if (atomlabels(i).eq."C2'") then
1803:                         dihed(2,3)=i1553:                         dihed(2,3)=i
1804:                                 dihed(3,2)=i1554:                                 dihed(3,2)=i
1805:         else if (atomlabels(i).eq."H2'1") then1555:         else if (atomlabels(i).eq."H2'1") then
1806: !                        dihed(2,4)=i            ! H3'-C3'-C2'-H2' should be cis1556: !                        dihed(2,4)=i            ! H3'-C3'-C2'-H2' should be cis
1807: !       !                         dihed(3,1)=i1557: !                                dihed(3,1)=i
1808:                 atom(5)=i1558:                 atom(5)=i
1809:         else if (atomlabels(i).eq."C1'") then1559:         else if (atomlabels(i).eq."C1'") then
1810:                                 dihed(3,3)=i1560:                                 dihed(3,3)=i
1811:         else if (atomlabels(i).eq."H1'") then1561:         else if (atomlabels(i).eq."H1'") then
1812:                                 dihed(3,4)=i    ! H2'-C2'-C1'-H1' should be trans1562:                                 dihed(3,4)=i    ! H2'-C2'-C1'-H1' should be trans
1813:                 atom(8)=i1563:                 atom(8)=i
1814:                 atom(7)=i+1 ! the N of the base follows the H1' atom in the topology file1564:                 atom(7)=i+1 ! the N of the base follows the H1' atom in the topology file
1815:         else if (atomlabels(i).eq."O2'") then1565:         else if (atomlabels(i).eq."O2'") then
1816:                 atom(6)=i1566:                 atom(6)=i
1817:         else if (atomlabels(i).eq."O3'") then1567:         else if (atomlabels(i).eq."O3'") then
1818:                 atom(4)=i1568:                 atom(4)=i
1819:                 resid=resid+11569:                 resid=resid+1
1820: !               do j=1,31570: !               do j=1,3
1821: !                call calc_dihedral(dihed(j,:),angle(j),coords,natoms) ! O3' is the last atom within one residue1571: !                call calc_dihedral(dihed(j,:),angle(j),coords,natoms) ! O3' is the last atom within one residue
1822: !               end do1572: !               end do
1823: 1573: 
1824:                call calc_distance(atom(1),atom(3),distance(1),coords,natoms) ! distances between C5' and H5' - H3' and O3'1574:                call calc_distance(atom(1),atom(3),distance(1),coords,natoms) ! distances between C5' and H5' - H3' and O3'
1825:                call calc_distance(atom(1),atom(4),distance(2),coords,natoms)1575:                call calc_distance(atom(1),atom(4),distance(2),coords,natoms)
1826:                call calc_distance(atom(2),atom(3),distance(4),coords,natoms)1576:                call calc_distance(atom(2),atom(3),distance(4),coords,natoms)
1827:                call calc_distance(atom(2),atom(4),distance(3),coords,natoms)1577:                call calc_distance(atom(2),atom(4),distance(3),coords,natoms)
1828: !              ! call calc_distance(atom(3),atom(5),distance(5),coords,natoms) ! not used, since this is DNA1578: !               call calc_distance(atom(3),atom(5),distance(5),coords,natoms) ! not used, since this is DNA
1829: !               call calc_distance(atom(3),atom(6),distance(6),coords,natoms)1579: !               call calc_distance(atom(3),atom(6),distance(6),coords,natoms)
1830: !               call calc_distance(atom(4),atom(5),distance(8),coords,natoms)1580: !               call calc_distance(atom(4),atom(5),distance(8),coords,natoms)
1831: !               call calc_distance(atom(4),atom(6),distance(7),coords,natoms)1581: !               call calc_distance(atom(4),atom(6),distance(7),coords,natoms)
1832: !               call calc_distance(atom(3),atom(7),distance(5),coords,natoms) ! distances between H1' and N1' - H3' and O3'1582: !               call calc_distance(atom(3),atom(7),distance(5),coords,natoms) ! distances between H1' and N1' - H3' and O3'
1833: !               call calc_distance(atom(3),atom(8),distance(6),coords,natoms)1583: !               call calc_distance(atom(3),atom(8),distance(6),coords,natoms)
1834:                call calc_distance(atom(4),atom(8),distance(5),coords,natoms)1584:                call calc_distance(atom(4),atom(8),distance(5),coords,natoms)
1835:                call calc_distance(atom(4),atom(7),distance(6),coords,natoms)1585:                call calc_distance(atom(4),atom(7),distance(6),coords,natoms)
1836:                call calc_distance(atom(2),atom(8),distance(7),coords,natoms)1586:                call calc_distance(atom(2),atom(8),distance(7),coords,natoms)
1837:                call calc_distance(atom(2),atom(9),distance(8),coords,natoms)1587:                call calc_distance(atom(2),atom(9),distance(8),coords,natoms)
1838:                call calc_distance(atom(2),atom(8),distance(9),coords,natoms)1588:                call calc_distance(atom(2),atom(8),distance(9),coords,natoms)
1839:                call calc_distance(atom(2),atom(10),distance(10),coords,natoms)1589:                call calc_distance(atom(2),atom(10),distance(10),coords,natoms)
1840:               !  write(*,*) distance(7:8)1590:                 write(*,*) distance(7:8)
1841:                 dihed(1,1)=i1591:                 dihed(1,1)=i
1842: 1592: 
1843:         do j=1,51593:         do j=1,5
1844:                 if(distance(2*j-1)>distance(2*j)+0.2D0) then1594:                 if(distance(2*j-1)>distance(2*j)+0.2D0) then
1845:                         goodstructure=.false.1595:                         goodstructure=.false.
1846:                    write(MYUNITNEW,'(A,I3,A,2I3,3F8.3)') 'residue ',resid,&1596:                    write(MYUNITNEW,'(A,I3,A,2I3,3F8.3)') 'residue ',resid,&
1847:                         &' wrong distances ', 2*j-1, 2*j,distance(2*j-1),distance(2*j)1597:                         &' wrong distances ', 2*j-1, 2*j,distance(2*j-1),distance(2*j)
1848: 1598: 
1849:                 end if1599:                 end if
1850:         end do1600:         end do
1851: !        if(abs(angle(1))<80.0D0 .or. abs(angle(2))>90.0D0 .or. abs(angle(3))<80.0D0) then1601: !        if(abs(angle(1))<80.0D0 .or. abs(angle(2))>90.0D0 .or. abs(angle(3))<80.0D0) then
1852: !            write(MYUNITNEW,'(A,I6,A,I3,A,3F8.3)') 'minimum ',k,' residue ',resid,' dihedrals: ',angle(1:3)1602: !            write(MYUNITNEW,'(A,I6,A,I3,A,3F8.3)') 'minimum ',k,' residue ',resid,' dihedrals: ',angle(1:3)
1853: !            goodstructure(k)=.false.1603: !            goodstructure(k)=.false.
1854: !        end if1604: !        end if
1855:        end if1605:        end if
1856:       end do1606:       end do
1857:         if (goodstructure) then1607:         if (goodstructure) then
1858:             write(MYUNITNEW,'(A)') 'Structure OK'1608:             write(MYUNITNEW,'(A)') 'structure OK'
1859:         end if1609:         end if
1860: 1610: 
1861: end subroutine check_cistrans_dna1611: end subroutine check_cistrans_dna
1862: 1612: 
1863: subroutine check_cistrans_protein(coords,natoms,goodstructure,minomega,cisarray)1613: subroutine check_cistrans_protein(coords,natoms,goodstructure,minomega,cisarray)
1864: use modamber9,only : myunitnew,i02,m02,m04,ih,ix1614: use modamber9,only : myunitnew,i02,m02,m04,ih,ix
1865: implicit none1615: implicit none
1866: 1616: 
1867: ! stuff for detecting cis-trans isomerisation of omega angles in peptide bonds 1617: ! stuff for detecting cis-trans isomerisation of omega angles in peptide bonds 
1868: logical,intent(out) :: goodstructure1618: logical,intent(out) :: goodstructure
1878: angle(:) = 0.0D01628: angle(:) = 0.0D0
1879: goodstructure=.true.1629: goodstructure=.true.
1880: cisarray(:)=01630: cisarray(:)=0
1881: DO i1=1,natoms1631: DO i1=1,natoms
1882: 1632: 
1883:  IF(ix(resnumber+i02-1)==i1) THEN1633:  IF(ix(resnumber+i02-1)==i1) THEN
1884:     currentresidue=resnumber1634:     currentresidue=resnumber
1885:     resnumber = resnumber + 11635:     resnumber = resnumber + 1
1886:  END IF1636:  END IF
1887: 1637: 
1888:  IF(ih(m04+i1-1)=='C   '.and.ih(m04+i1)=='O   '.and.ih(m04+i1gv+1)=='N   '.and.ih(m04+i1+2)=='H   ') THEN1638:  IF(ih(m04+i1-1)=='C   '.and.ih(m04+i1)=='O   '.and.ih(m04+i1+1)=='N   '.and.ih(m04+i1+2)=='H   ') THEN
1889:       dihed(1,2)=i11639:       dihed(1,2)=i1
1890:  ELSE IF(ih(m04+i1-1)=='O   '.and.ih(m04+i1-2)=='C   '.and.ih(m04+i1)=='N   '.and.ih(m04+i1+1)=='H   ') THEN1640:  ELSE IF(ih(m04+i1-1)=='O   '.and.ih(m04+i1-2)=='C   '.and.ih(m04+i1)=='N   '.and.ih(m04+i1+1)=='H   ') THEN
1891:       dihed(1,1)=i11641:       dihed(1,1)=i1
1892:  ELSE IF     (ih(m04+i1-1)=='N   '.and.ih(m04+i1-2)=='O   '.and.ih(m04+i1-3)=='C   '.and.ih(m04+i1)=='H   ') THEN1642:  ELSE IF     (ih(m04+i1-1)=='N   '.and.ih(m04+i1-2)=='O   '.and.ih(m04+i1-3)=='C   '.and.ih(m04+i1)=='H   ') THEN
1893:       dihed(1,3)=i11643:       dihed(1,3)=i1
1894:  ELSE IF((ih(m04+i1-1)=='H   '.OR.ih(m04+i1-1)=='CD  ').AND.&1644:  ELSE IF((ih(m04+i1-1)=='H   '.OR.ih(m04+i1-1)=='CD  ').AND.&
1895:         & ih(m04+i1-2)=='N   '.AND.ih(m04+i1-3)=='O   '.AND.ih(m04+i1-4)=='C   ') THEN 1645:         & ih(m04+i1-2)=='N   '.AND.ih(m04+i1-3)=='O   '.AND.ih(m04+i1-4)=='C   ') THEN 
1896: ! this is the last atom in one peptide bond in the topology file, so now1646: ! this is the last atom in one peptide bond in the topology file, so now
1897:                                     ! we can calculate the dihedral1647:                                     ! we can calculate the dihedral
1898:         dihed(1,4)=i11648:         dihed(1,4)=i1


r28580/egb.f 2017-03-30 13:30:32.588546696 +0100 r28579/egb.f 2017-03-30 13:30:34.400570786 +0100
704: !        for Ri=Rborn (assuming that rgbmax<cut)704: !        for Ri=Rborn (assuming that rgbmax<cut)
705: !   double precision spla,splb,splc,spld705: !   double precision spla,splb,splc,spld
706:    double precision rcut,rcutminusroff, rcutminusroff3inv,rcutminusroffinv!,dfgbdr706:    double precision rcut,rcutminusroff, rcutminusroff3inv,rcutminusroffinv!,dfgbdr
707: !   integer icount2707: !   integer icount2
708: 708: 
709:    common/extrasamberint/ifswitch,irespa2,nrespa2709:    common/extrasamberint/ifswitch,irespa2,nrespa2
710:    common/extrasamberdp/fswitchbeta710:    common/extrasamberdp/fswitchbeta
711: 711: 
712: !   logical calledonce712: !   logical calledonce
713: !   common/logicals/ calledonce713: !   common/logicals/ calledonce
714:  
715: !NaN test 
716: !OPEN(UNIT=234,FILE="NaN_debug") 
717:  
718:  
719:  
720: !============ QMMM ==========================================================714: !============ QMMM ==========================================================
721: !  If ifqnt == True and igb /=0 and igb /= 6 then we will do EGB with QMMM.715: !  If ifqnt == True and igb /=0 and igb /= 6 then we will do EGB with QMMM.
722: 716: 
723: !  In this case the charge array should currently be zero for the QM atoms.717: !  In this case the charge array should currently be zero for the QM atoms.
724: !  This corresponds to qmgb = 0 where we do only EGB(MM-MM) -- skipping 718: !  This corresponds to qmgb = 0 where we do only EGB(MM-MM) -- skipping 
725: !  QM-MM and QM-QM interactions.719: !  QM-MM and QM-QM interactions.
726: 720: 
727: !  If qmgb==1 then we need to copy back the qm resp charges from 721: !  If qmgb==1 then we need to copy back the qm resp charges from 
728: !  qm_resp_charges, so we can do EGB on everything. When done we need to 722: !  qm_resp_charges, so we can do EGB on everything. When done we need to 
729: !  ensure we zero the QM charge array again.723: !  ensure we zero the QM charge array again.
948:       icount = 0942:       icount = 0
949: !      icount2 = 0943: !      icount2 = 0
950: term2(1:natom) = one 944: term2(1:natom) = one 
951: term3(1:natom) = zero945: term3(1:natom) = zero
952: 946: 
953:       do j=i+1,natom947:       do j=i+1,natom
954:          xij = xi - x(3*j-2)948:          xij = xi - x(3*j-2)
955:          yij = yi - x(3*j-1)949:          yij = yi - x(3*j-1)
956:          zij = zi - x(3*j  )950:          zij = zi - x(3*j  )
957:          r2 = xij*xij + yij*yij + zij*zij951:          r2 = xij*xij + yij*yij + zij*zij
958: !dc550 To check distances952: 
959: !IF(r2.le.0.5)THEN953:          if( r2 <= cut .and. (onstep .or. r2 <= cut_inner)) then
960:  !WRITE(234,*)'The wrong distances are', 'ATOM',i,'ATOM',j, 'DISTANCE',r2 
961:   !ENDIF         
962:  if( r2 <= cut .and. (onstep .or. r2 <= cut_inner)) then 
963: 954: 
964: ! set for non-LES case so we can refer to reff_j not reff(j)955: ! set for non-LES case so we can refer to reff_j not reff(j)
965: ! this is needed with LES since each atom j has multiple reff956: ! this is needed with LES since each atom j has multiple reff
966: 957: 
967:              reff_j = reff(j)958:              reff_j = reff(j)
968:              icount = icount + 1959:              icount = icount + 1
969:              temp_jj(icount) = j960:              temp_jj(icount) = j
970:              r2x(icount) = r2 961:              r2x(icount) = r2 
971: 962: 
972: ! calculate term2 and term3 for force shifting963: ! calculate term2 and term3 for force shifting
2089:        f(3*mm_no  ) = f(3*mm_no  ) + dumz2080:        f(3*mm_no  ) = f(3*mm_no  ) + dumz
2090:      end do ! i = 1, qmmm_struct%nlink2081:      end do ! i = 1, qmmm_struct%nlink
2091:      call timer_stop(125)2082:      call timer_stop(125)
2092:      call timer_stop(103)2083:      call timer_stop(103)
2093: !------------ END MM LINK PAIR VDW ---------------------2084: !------------ END MM LINK PAIR VDW ---------------------
2094:      call timer_stop(23)2085:      call timer_stop(23)
2095:    end if !if ifqnt.2086:    end if !if ifqnt.
2096: !======== END QMMM ==========2087: !======== END QMMM ==========
2097: 2088: 
2098:    return2089:    return
2099:  
2100: !CLOSE(UNIT=234) 
2101: end subroutine egb 2090: end subroutine egb 
2102: 2091: 
2103: subroutine egb_calc_radii(igb,natom,x,fs,reff,onereff,fsmax,rgbmax, &2092: subroutine egb_calc_radii(igb,natom,x,fs,reff,onereff,fsmax,rgbmax, &
2104:                           rborn, offset,gbalpha,gbbeta,gbgamma,rbornstat, &2093:                           rborn, offset,gbalpha,gbbeta,gbgamma,rbornstat, &
2105:                           rbave,rbfluct,rbmax,rbmin, gbneckscale, ncopy, rdt &2094:                           rbave,rbfluct,rbmax,rbmin, gbneckscale, ncopy, rdt &
2106: 2095: 
2107: 2096: 
2108: 2097: 
2109:                           )2098:                           )
2110: 2099: 


r28580/ene.f 2017-03-30 13:30:32.868550421 +0100 r28579/ene.f 2017-03-30 13:30:35.756588783 +0100
317: 317: 
318: 318: 
319:    319:    
320:    !------------- local variables  ---------------------------320:    !------------- local variables  ---------------------------
321:    integer max190321:    integer max190
322:    parameter (max190=190)322:    parameter (max190=190)
323:    integer nb,ist,maxlen,ksum,i3,j3,ic,jn,ii,jj323:    integer nb,ist,maxlen,ksum,i3,j3,ic,jn,ii,jj
324:    double precision xij(max190),yij(max190),zij(max190),eaw(max190), &324:    double precision xij(max190),yij(max190),zij(max190),eaw(max190), &
325:          rij(max190),dfw(max190)325:          rij(max190),dfw(max190)
326:    double precision da,df,xa,ya,za,rij0,ebl326:    double precision da,df,xa,ya,za,rij0,ebl
327:   
328:   
329:    logical skip327:    logical skip
330:    integer piece,start,end,newnb328:    integer piece,start,end,newnb
331: 329: 
332: !! ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''330: !! ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
333: !! ;  EVB variables                                                   ;331: !! ;  EVB variables                                                   ;
334: !! ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,332: !! ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
335: 333: 
336:    integer :: ndx334:    integer :: ndx
337: 335: 
338: 336: 
377:          i3 = ib(jn+ist)375:          i3 = ib(jn+ist)
378:          j3 = jb(jn+ist)376:          j3 = jb(jn+ist)
379:          377:          
380:          !           ----- CALCULATION OF THE bond vector -----378:          !           ----- CALCULATION OF THE bond vector -----
381:        379:        
382:          xij(jn) = x(i3+1)-x(j3+1)380:          xij(jn) = x(i3+1)-x(j3+1)
383:          yij(jn) = x(i3+2)-x(j3+2)381:          yij(jn) = x(i3+2)-x(j3+2)
384:          zij(jn) = x(i3+3)-x(j3+3)382:          zij(jn) = x(i3+3)-x(j3+3)
385:       end do383:       end do
386:       384:       
387: !----------CHECK IF TWO ATOMS ARE GETTING TOO CLOSE, MIGHT CAUSE PROBLEMS IN INTERPOLATION--------385:       do jn = 1,maxlen
388: !dc550 
389:  !  do jn = 1,maxlen 
390:   !    i3 = ib(jn+ist) 
391:     !  j3 = jb(jn+ist) 
392:   !    rij0 = xij(jn)*xij(jn)+yij(jn)*yij(jn)+zij(jn)*zij(jn) 
393:   !    rij(jn) = sqrt(rij0) 
394:   !    IF(rij(jn).le.0.5)THEN 
395:   !    xij(jn) = 0. 
396:   !    yij(jn) = 0. 
397:   !    zij(jn) = 0. 
398:   !    x(i3+1) = x(i3+1)+0.3 
399:   !    x(j3+1) = x(j3+1)-0.3 
400:   !    x(i3+2) = x(i3+2)+0.3 
401:   !    x(j2+2) = x(j3+2)-0.3 
402:   !    x(i3+3) = x(i3+3)+0.3 
403: ! 
404:  !     yij(jn) = x(i3+2) - x(j3+2) 
405:   !    zij(jn) = x(i3+3) - x(j3+3) 
406:   !    !swap values 
407:   !    ENDIF        
408:  ! enddo 
409:  
410:    do jn = 1,maxlen 
411:          rij0 = xij(jn)*xij(jn)+yij(jn)*yij(jn)+zij(jn)*zij(jn)386:          rij0 = xij(jn)*xij(jn)+yij(jn)*yij(jn)+zij(jn)*zij(jn)
412:          rij(jn) = sqrt(rij0)387:          rij(jn) = sqrt(rij0)
413: 388:       end do
414:     end do 
415:       !         ----- CALCULATION OF THE ENERGY AND DER -----389:       !         ----- CALCULATION OF THE ENERGY AND DER -----
416:       390:       
417:       do jn = 1,maxlen391:       do jn = 1,maxlen
418:          ic = icb(jn+ist)392:          ic = icb(jn+ist)
419:          rij0 = rij(jn)393:          rij0 = rij(jn)
420:          da = rij0-req(ic)394:          da = rij0-req(ic)
421:          !                                 for rms deviation from ideal bonds:395:          !                                 for rms deviation from ideal bonds:
422:          ebdev = ebdev + da*da396:          ebdev = ebdev + da*da
423:          df = rk(ic)*da397:          df = rk(ic)*da
424:          eaw(jn) = df*da398:          eaw(jn) = df*da


r28580/nonbon.f 2017-03-30 13:30:33.212555158 +0100 r28579/nonbon.f 2017-03-30 13:30:44.884710098 +0100
 52:          i3 = 3*i-3 52:          i3 = 3*i-3
 53:          cgi = cg(i) 53:          cgi = cg(i)
 54:          iaci = iac(i) 54:          iaci = iac(i)
 55:          jaci = ntypes*(iaci-1) 55:          jaci = ntypes*(iaci-1)
 56:          do 50 jj = 1,npr 56:          do 50 jj = 1,npr
 57:             lpair = lpair+1 57:             lpair = lpair+1
 58:             j = iar2(lpair) 58:             j = iar2(lpair)
 59:             j3 = 3*j-3 59:             j3 = 3*j-3
 60:             do 10 m = 1,3 60:             do 10 m = 1,3
 61:                xij(m) = x(i3+m)-x(j3+m) 61:                xij(m) = x(i3+m)-x(j3+m)
 62:               62:             10 continue
 63:               63:             s = xij(1)**2 + xij(2)**2 + xij(3)**2
 64:               64:             if (ntb /= 0) then
 65:   
 66:            10 continue 
 67: !dc550 addition 
 68:                      
 69:        !IF(ntb==0)THEN 
 70:        !  DO m = 1,3 
 71:         !  IF(s.le.0.5)THEN 
 72:          !     x(i3+m) = x(i3+m)+0.3 
 73:          !     x(j3+m) = x(j3+m)-0.3 
 74:          ! ENDIF 
 75:           ! ENDDO 
 76:        ! ENDIF 
 77:        s = xij(1)**2+xij(2)**2+xij(3)**2 
 78:   
 79:    if (ntb /= 0) then 
 80:                do 20 m = 1,3 65:                do 20 m = 1,3
 81:                   if(xij(m) >= boxh(m)) then 66:                   if(xij(m) >= boxh(m)) then
 82:                      xij(m) = xij(m)-box(m) 67:                      xij(m) = xij(m)-box(m)
 83:                   else if(xij(m) < -boxh(m)) then 68:                   else if(xij(m) < -boxh(m)) then
 84:                      xij(m) = xij(m)+box(m) 69:                      xij(m) = xij(m)+box(m)
 85:                   end if 70:                   end if
 86:                20 continue 71:                20 continue
 87:                s = xij(1)**2+xij(2)**2+xij(3)**2 72:                s = xij(1)**2+xij(2)**2+xij(3)**2
 88:                if(ntb < 0) then 73:                if(ntb < 0) then
 89:                   s = s+dmin1(0.0d0,boxoq-dabs(xij(1))-dabs(xij(2)) & 74:                   s = s+dmin1(0.0d0,boxoq-dabs(xij(1))-dabs(xij(2)) &
 90:                         -dabs(xij(3)))*box(1) 75:                         -dabs(xij(3)))*box(1)
 91:                else if(ntm /= 0) then 76:                else if(ntm /= 0) then
 92:         77:                   s = s+cosb2*xij(1)*xij(3)
 93:            s = s+cosb2*xij(1)*xij(3) 
 94:                end if 78:                end if
 95:             end if 79:             end if
 96:              80:             
 97:             iacj = iac(j) 81:             iacj = iac(j)
 98:             index = jaci+iacj 82:             index = jaci+iacj
 99:             ic = ico(index) 83:             ic = ico(index)
100:             if(ntm /= 0) call traco(1,0,xij,beta,-1) 84:             if(ntm /= 0) call traco(1,0,xij,beta,-1)
101:              85:             
102:             !     ----- claculate the electrostaic energy ----- 86:             !     ----- claculate the electrostaic energy -----
103:              87:             


r28580/printe.f 2017-03-30 13:30:33.500558819 +0100 r28579/printe.f 2017-03-30 13:30:45.248714952 +0100
858:    gradient_max = ZERO858:    gradient_max = ZERO
859:    atom_number_of_gmax = 1859:    atom_number_of_gmax = 1
860:    do i = 1,3*natom860:    do i = 1,3*natom
861:       gi = abs(gradient(i))861:       gi = abs(gradient(i))
862:       if (gi > gradient_max) then862:       if (gi > gradient_max) then
863:          gradient_max = gi863:          gradient_max = gi
864:          atom_number_of_gmax = i864:          atom_number_of_gmax = i
865:       end if865:       end if
866:    end do866:    end do
867:    atom_number_of_gmax = (atom_number_of_gmax - 1)/3 + 1867:    atom_number_of_gmax = (atom_number_of_gmax - 1)/3 + 1
868: !> DUMP THE GRADIENTS HERE 
869:    OPEN(UNIT=109, FILE="amber_grad", STATUS="unknown") 
870:  
871:    WRITE(109,*) "BEFORE STEP******* CHECKPOINT1" 
872:    DO i = 1,3*natom 
873:    WRITE(109,*) gradient(i) 
874:    ENDDO 
875:    WRITE(109,*) "END STEP*********CHECKPOINT1" 
876:    CLOSE(109) 
877:  
878:  
879:  
880:  
881:  
882: 868: 
883:    return869:    return
884: end subroutine grdmax870: end subroutine grdmax
885: 871: 
886: 872: 
887: !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++873: !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
888: !+ Print status of a minimization calculation.874: !+ Print status of a minimization calculation.
889: !-----------------------------------------------------------------------875: !-----------------------------------------------------------------------
890: !     --- PRINTE ---876: !     --- PRINTE ---
891: !-----------------------------------------------------------------------877: !-----------------------------------------------------------------------
1114:    dvdl   = energies(17)1100:    dvdl   = energies(17)
1115:    diprms = energies(22)1101:    diprms = energies(22)
1116:    dipiter = energies(23)1102:    dipiter = energies(23)
1117:    dipole_temp = energies(24)1103:    dipole_temp = energies(24)
1118:    escf = energies(25)1104:    escf = energies(25)
1119:    edisp  = energies(28)1105:    edisp  = energies(28)
1120: 1106: 
1121: !   write(6,9018)1107: !   write(6,9018)
1122:    write(6,9028) nstep, epot, gradient_rms, gradient_max, &1108:    write(6,9028) nstep, epot, gradient_rms, gradient_max, &
1123:          atom_name_of_gmax, atom_number_of_gmax1109:          atom_name_of_gmax, atom_number_of_gmax
1124: !  write(6,9028)  
1125: ! sf344> commented out1110: ! sf344> commented out
1126: !   write(6,9038) ebond,eangle,edihed1111: !   write(6,9038) ebond,eangle,edihed
1127:    if( igb == 0 ) then1112:    if( igb == 0 ) then
1128: !      write(6,9048) enonb,enele,ehbond1113: !      write(6,9048) enonb,enele,ehbond
1129:    else if ( igb == 10 ) then1114:    else if ( igb == 10 ) then
1130: !      write(6,9050) enonb,enele,ehbond1115: !      write(6,9050) enonb,enele,ehbond
1131:    else1116:    else
1132: !      write(6,9049) enonb,enele,ehbond1117: !      write(6,9049) enonb,enele,ehbond
1133:    end if1118:    end if
1134: !   write(6,9058) enb14,eel14,econst1119: !   write(6,9058) enb14,eel14,econst


r28580/veclib.f 2017-03-30 13:30:33.784562598 +0100 r28579/veclib.f 2017-03-30 13:30:45.524718626 +0100
115: end subroutine vdln 115: end subroutine vdln 
116: !--------------------------------------------------------------116: !--------------------------------------------------------------
117: 117: 
118: !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++118: !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
119: !+ vectorized inverse square root119: !+ vectorized inverse square root
120: subroutine vdinvsqrt( n, x, y )120: subroutine vdinvsqrt( n, x, y )
121: 121: 
122:    use constants, only : one 122:    use constants, only : one 
123: 123: 
124:    implicit none124:    implicit none
125:    integer  n, i125:    integer  n
126:    double precision   x(n), y(n)126:    double precision   x(n), y(n)
127: !   OPEN(UNIT=110,FILE="dump_sqrt")  127:    
128:  128: 
129:  !     WRITE(110,*) "--------------------------------------"129: 
130:   !    WRITE(110,*) (sqrt(x(i)),i=1,n)130: 
131:       y(1:n) = one/sqrt( x(1:n) )131:       y(1:n) = one/sqrt( x(1:n) )
132: 132: 
133:    133:    
134:    return134:    return
135:    !CLOSE(UNIT=110) 
136: end subroutine vdinvsqrt 135: end subroutine vdinvsqrt 
137: !--------------------------------------------------------------136: !--------------------------------------------------------------
138: 137: 
139: !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++138: !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
140: !+ vectorized inverse139: !+ vectorized inverse
141: subroutine vdinv( n, x, y )140: subroutine vdinv( n, x, y )
142: 141: 
143:    use constants, only : one 142:    use constants, only : one 
144: 143: 
145:    implicit none144:    implicit none


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0