hdiff output

r22833/amberinterface.f 2017-03-30 13:30:54.128832968 +0100 r22832/amberinterface.f 2017-03-30 13:30:54.416836794 +0100
667: use nblist, only : nbflag, skinnb667: use nblist, only : nbflag, skinnb
668: use porfuncs668: use porfuncs
669: 669: 
670: implicit none670: implicit none
671: 671: 
672: 672: 
673: integer JP, itime1, now(3), nmovableatoms,i,j,movableatomlist(nmovableatoms),nblocks,atomsinblock(nblocks),offset1,k673: integer JP, itime1, now(3), nmovableatoms,i,j,movableatomlist(nmovableatoms),nblocks,atomsinblock(nblocks),offset1,k
674: double precision        :: y(3*natoms), grad(3*natoms), ereal,ligandcentre(3),ligandcoords(3),ligandcoordsrotated(3)674: double precision        :: y(3*natoms), grad(3*natoms), ereal,ligandcentre(3),ligandcoords(3),ligandcoordsrotated(3)
675: double precision        :: rotationmatrix(nblocks,3,3),twopi,pi,DPRAND,randomphi,randompsi,randomtheta,DIST, dummyz675: double precision        :: rotationmatrix(nblocks,3,3),twopi,pi,DPRAND,randomphi,randompsi,randomtheta,DIST, dummyz
676: double precision        :: st,ct,sph,cph,sps,cps, VECBAX, VECBAY, VECBAZ, randomx, randomy, randomz, dummyx, dummyy676: double precision        :: st,ct,sph,cph,sps,cps, VECBAX, VECBAY, VECBAZ, randomx, randomy, randomz, dummyx, dummyy
677: double precision        :: atomblockcentre(nblocks,3),distancematrix(natoms,natoms)677: double precision        :: atomblockcentre(nblocks,3)
678: logical                 :: ligmovet,mdstept1,randomseedt2,reschangedt(nres),blockmovet,overlapt678: logical                 :: ligmovet,mdstept1,randomseedt2,reschangedt(nres),blockmovet
679: character(len=10)       :: datechar,timechar,zonechar679: character(len=10)       :: datechar,timechar,zonechar
680: integer                 :: values(8),iostatus,restemp,i1,currentresidue,resnumber,overlaparray(natoms,natoms)680: integer                 :: values(8),iostatus,restemp,i1,currentresidue,resnumber
681: integer                 :: overlaparraysave(natoms,natoms) 
682: character(len=10)       :: rotmaxchangestr,rotpselectstr,rotcutoffstr,rotcentrestr,rotoccuwstr681: character(len=10)       :: rotmaxchangestr,rotpselectstr,rotcutoffstr,rotcentrestr,rotoccuwstr
683: ! parameters682: ! parameters
684: pi=ATAN(1.0D0)*4683: pi=ATAN(1.0D0)*4
685: twopi=2.0D0*pi684: twopi=2.0D0*pi
686: ! increments with each takestepamber call685: ! increments with each takestepamber call
687: n_amb_calls=n_amb_calls+1686: n_amb_calls=n_amb_calls+1
688: if(nmovableatoms==0.and.(ligmovet.or.blockmovet)) then 687: if(nmovableatoms==0.and.(ligmovet.or.blockmovet)) then 
689:         write(MYUNITNEW,*)&688:         write(MYUNITNEW,*)&
690:         &        ' takestepamber> ERROR - need to specify movable atoms and use the MOVABLEATOMS keyword!'689:         &        ' takestepamber> ERROR - need to specify movable atoms and use the MOVABLEATOMS keyword!'
691:         STOP690:         STOP
743: 742: 
744: ! LIGAND MOVES START HERE 743: ! LIGAND MOVES START HERE 
745: ! At the moment we have rigid rotation and cartesian perturbation. There is also a translation for use in744: ! At the moment we have rigid rotation and cartesian perturbation. There is also a translation for use in
746: ! steered minimisation745: ! steered minimisation
747: 746: 
748: ! 1) RIGID BODY ROTATION747: ! 1) RIGID BODY ROTATION
749: 748: 
750: IF((LIGMOVET.OR.STEEREDMINT).AND.doligmove) THEN749: IF((LIGMOVET.OR.STEEREDMINT).AND.doligmove) THEN
751: 750: 
752:         IF (STEEREDMINT) WRITE(MYUNITNEW,'(A)') 'takestepamber> rotating and translating the ligand'751:         IF (STEEREDMINT) WRITE(MYUNITNEW,'(A)') 'takestepamber> rotating and translating the ligand'
753:   overlapt=.true. 
754:   overlaparraysave(:,:)=0 
755:   overlaparray(:,:)=0 
756:   distancematrix(:,:)=0.0D0 
757:         IF(BLOCKMOVET) THEN 
758:          ! Try to avoid moves that cause extremely high initial forces (or cold fusion), 
759:          ! by redoing the rotation moves whenever two atoms get closer than 0.5 Angstroms. 
760:           do i=1,natoms 
761:            do j=1,natoms 
762:               if(i==j) cycle 
763:               distancematrix(i,j)=sqrt((y(3*i-2)-y(3*j-2))**2+(y(3*i-1)-y(3*j-1))**2+(y(3*i)-y(3*j))**2) 
764:               if(distancematrix(i,j)<0.7D0) then  
765:                 overlaparraysave(i,j)=1  !there's an initial overlap 
766: !                write(*,*) 'atoms close together before rotating ', i,j,distancematrix(i,j) 
767:               end if 
768:            end do 
769:           end do 
770:         END IF 
771:  DO WHILE(OVERLAPT) 
772:      overlapt=.false. 
773:      DO i=1,nblocks ! nblocks=1 if BLOCKMOVET=.FALSE. 752:      DO i=1,nblocks ! nblocks=1 if BLOCKMOVET=.FALSE. 
774:         randomphi=(DPRAND()-0.5)*twopi*ligrotscale753:         randomphi=(DPRAND()-0.5)*twopi*ligrotscale
775:         randomtheta=(DPRAND()-0.5)*pi*ligrotscale754:         randomtheta=(DPRAND()-0.5)*pi*ligrotscale
776:         randompsi=(DPRAND()-0.5)*twopi*ligrotscale755:         randompsi=(DPRAND()-0.5)*twopi*ligrotscale
777:         st=sin(randomtheta)756:         st=sin(randomtheta)
778:         ct=cos(randomtheta)757:         ct=cos(randomtheta)
779:         sph=sin(randomphi)758:         sph=sin(randomphi)
780:         cph=cos(randomphi)759:         cph=cos(randomphi)
781:         sps=sin(randompsi)760:         sps=sin(randompsi)
782:         cps=cos(randompsi)761:         cps=cos(randompsi)
795:      END DO774:      END DO
796: !         y(1:3*natom)=coords1(1:3*natom)775: !         y(1:3*natom)=coords1(1:3*natom)
797: 776: 
798: !        y(1:3*natom)=x(lcrd:lcrd+3*natom-1)777: !        y(1:3*natom)=x(lcrd:lcrd+3*natom-1)
799: 778: 
800: ! work out the centre of coordinates for the ligand779: ! work out the centre of coordinates for the ligand
801: 780: 
802:         IF(BLOCKMOVET) THEN781:         IF(BLOCKMOVET) THEN
803:           atomblockcentre(:,:)=0.0D0782:           atomblockcentre(:,:)=0.0D0
804:           offset1=0783:           offset1=0
805: !          overlapt=.false. 
806:           do i=1,nblocks784:           do i=1,nblocks
807:              if(i>1) offset1=offset1+atomsinblock(i-1)785:              if(i>1) offset1=offset1+atomsinblock(i-1)
808:              do k=1,atomsinblock(i)786:              do k=1,atomsinblock(i)
809:                 j=movableatomlist(k+offset1)787:                 j=movableatomlist(k+offset1)
810: !                WRITE(*,*) k, j788: !                WRITE(*,*) k, j
811:                 atomblockcentre(i,1)=atomblockcentre(i,1)+y(3*j-2)789:                 atomblockcentre(i,1)=atomblockcentre(i,1)+y(3*j-2)
812:                 atomblockcentre(i,2)=atomblockcentre(i,2)+y(3*j-1)790:                 atomblockcentre(i,2)=atomblockcentre(i,2)+y(3*j-1)
813:                 atomblockcentre(i,3)=atomblockcentre(i,3)+y(3*j  )791:                 atomblockcentre(i,3)=atomblockcentre(i,3)+y(3*j  )
814:              end do792:              end do
815:                 atomblockcentre(i,:)=atomblockcentre(i,:)/atomsinblock(i)793:                 atomblockcentre(i,:)=atomblockcentre(i,:)/atomsinblock(i)
820:                 ligandcoords(2)=y(3*j-1)-atomblockcentre(i,2)798:                 ligandcoords(2)=y(3*j-1)-atomblockcentre(i,2)
821:                 ligandcoords(3)=y(3*j  )-atomblockcentre(i,3)799:                 ligandcoords(3)=y(3*j  )-atomblockcentre(i,3)
822: ! rotate the block of atoms with random rotation matrix800: ! rotate the block of atoms with random rotation matrix
823:                 ligandcoordsrotated=MATMUL(rotationmatrix(i,:,:),ligandcoords)801:                 ligandcoordsrotated=MATMUL(rotationmatrix(i,:,:),ligandcoords)
824: ! translate back the block of atoms to the original centre of their coordinates802: ! translate back the block of atoms to the original centre of their coordinates
825:                 y(3*j-2)=ligandcoordsrotated(1)+atomblockcentre(i,1)803:                 y(3*j-2)=ligandcoordsrotated(1)+atomblockcentre(i,1)
826:                 y(3*j-1)=ligandcoordsrotated(2)+atomblockcentre(i,2)804:                 y(3*j-1)=ligandcoordsrotated(2)+atomblockcentre(i,2)
827:                 y(3*j  )=ligandcoordsrotated(3)+atomblockcentre(i,3)805:                 y(3*j  )=ligandcoordsrotated(3)+atomblockcentre(i,3)
828:              end do            806:              end do            
829:           end do807:           end do
830:           ! compute the new distances between all atoms after the rotation move and 
831:           ! determine whether the overlap array has changed. 
832:           do i=1,natoms 
833:            do j=1,natoms 
834:               if(i==j) cycle 
835:               distancematrix(i,j)=sqrt((y(3*i-2)-y(3*j-2))**2+(y(3*i-1)-y(3*j-1))**2+(y(3*i)-y(3*j))**2) 
836:               if(distancematrix(i,j)<0.7D0)then 
837:                  overlaparray(i,j)=1 
838:               else 
839:                  overlaparray(i,j)=0  
840:               end if 
841:               if(overlaparraysave(i,j)-overlaparray(i,j)<0) then  
842:                 overlapt=.true. 
843: !                write(*,*) 'atoms close together after rotating ', i,j,distancematrix(i,j) 
844:                 GOTO 95 
845:               end if 
846:            end do 
847:          end do 
848:         ELSE808:         ELSE
849:           ligandcentre(:)=0.0D0809:           ligandcentre(:)=0.0D0
850:           do i=1,nmovableatoms810:           do i=1,nmovableatoms
851:                 j=movableatomlist(i)811:                 j=movableatomlist(i)
852:                 ligandcentre(1)=ligandcentre(1)+y(3*j-2)812:                 ligandcentre(1)=ligandcentre(1)+y(3*j-2)
853:                 ligandcentre(2)=ligandcentre(2)+y(3*j-1)813:                 ligandcentre(2)=ligandcentre(2)+y(3*j-1)
854:                 ligandcentre(3)=ligandcentre(3)+y(3*j  )814:                 ligandcentre(3)=ligandcentre(3)+y(3*j  )
855:           end do815:           end do
856:                 ligandcentre(:) = ligandcentre(:)/nmovableatoms816:                 ligandcentre(:) = ligandcentre(:)/nmovableatoms
857: !        write(MYUNITNEW,*) 'ligandcentre = ', ligandcentre(:)817: !        write(MYUNITNEW,*) 'ligandcentre = ', ligandcentre(:)
865:                 ligandcoordsrotated=MATMUL(rotationmatrix(1,:,:),ligandcoords)825:                 ligandcoordsrotated=MATMUL(rotationmatrix(1,:,:),ligandcoords)
866: !          write(MYUNITNEW,*) 'rotated coords: ', ligandcoordsrotated(:)826: !          write(MYUNITNEW,*) 'rotated coords: ', ligandcoordsrotated(:)
867: !          write(MYUNITNEW,*) 'original coords: ', ligandcoords(:)827: !          write(MYUNITNEW,*) 'original coords: ', ligandcoords(:)
868: 828: 
869: ! translate the ligand atoms back to the centre of the original ligand                829: ! translate the ligand atoms back to the centre of the original ligand                
870:                 y(3*j-2)=ligandcoordsrotated(1)+ligandcentre(1)830:                 y(3*j-2)=ligandcoordsrotated(1)+ligandcentre(1)
871:                 y(3*j-1)=ligandcoordsrotated(2)+ligandcentre(2)831:                 y(3*j-1)=ligandcoordsrotated(2)+ligandcentre(2)
872:                 y(3*j  )=ligandcoordsrotated(3)+ligandcentre(3)832:                 y(3*j  )=ligandcoordsrotated(3)+ligandcentre(3)
873:           end do833:           end do
874:         END IF834:         END IF
875: 95        CONTINUE835: 
876:  END DO  ! while(overlapt) 
877: ! 2) RANDOM CARTESIAN PERTURBATION836: ! 2) RANDOM CARTESIAN PERTURBATION
878:         if ((LIGMOVET).and.(ligcartstep.gt.0)) then837:         if ((LIGMOVET).and.(ligcartstep.gt.0)) then
879:            do i=1,nmovableatoms838:            do i=1,nmovableatoms
880:               j=movableatomlist(i)839:               j=movableatomlist(i)
881:               randomx=DPRAND()840:               randomx=DPRAND()
882:               randomy=DPRAND()841:               randomy=DPRAND()
883:               randomz=DPRAND()842:               randomz=DPRAND()
884:               y(3*j-2)=y(3*j-2)+2*randomx*ligcartstep-ligcartstep843:               y(3*j-2)=y(3*j-2)+2*randomx*ligcartstep-ligcartstep
885:               y(3*j-1)=y(3*j-1)+2*randomy*ligcartstep-ligcartstep844:               y(3*j-1)=y(3*j-1)+2*randomy*ligcartstep-ligcartstep
886:               y(3*j  )=y(3*j  )+2*randomz*ligcartstep-ligcartstep845:               y(3*j  )=y(3*j  )+2*randomz*ligcartstep-ligcartstep
1635:  & WARNING -- cis (or deformed) peptide bond detected between residues ',&1594:  & WARNING -- cis (or deformed) peptide bond detected between residues ',&
1636:  & currentresidue-1, currentresidue1595:  & currentresidue-1, currentresidue
1637:         END IF1596:         END IF
1638:       END IF1597:       END IF
1639: !         WRITE(*,*) currentresidue-1,ih(currentresidue-1+m02-1), angle(1), minomega1598: !         WRITE(*,*) currentresidue-1,ih(currentresidue-1+m02-1), angle(1), minomega
1640:  END IF1599:  END IF
1641: 1600: 
1642: END DO1601: END DO
1643: end subroutine check_cistrans_protein1602: end subroutine check_cistrans_protein
1644: 1603: 
1645: subroutine check_chirality_generic(coords,natoms,goodstructure,init) 
1646: use modamber9,only : myunitnew,i02,m02,m04,ih,ix,uachiral,dihedralsave,atomindex,exclude 
1647: implicit none 
1648: ! Checks the chirality of each carbon atom, and compares it to that of the initial structure. 
1649: ! Returns goodstructure=.false. if there is a change in sign of the calculated dihedrals. 
1650: ! Does not check chirality if two hydrogens are attached to the same carbon. 
1651: ! Works by finding the four closest atoms to the carbon atom in question, ordering the  
1652: ! indices of those ascending, and calculating the dihedral. 
1653:  
1654: integer,intent(in) :: natoms 
1655: double precision,intent(in) :: coords(3*natoms) 
1656: integer :: i,i1,i2,i3,j,currentresidue,resnumber,hcount 
1657: double precision :: dihedral(natoms),distancematrix(natoms,natoms),dist 
1658: logical :: goodstructure,init,sp(natoms),sp2(natoms) 
1659: integer :: chiarray(natoms) 
1660: character(len=1) :: labelchar(natoms) 
1661:  
1662: goodstructure=.true. 
1663:  
1664: ! WRITE(*,*) 'init=', init 
1665: if(init) then 
1666:  exclude(:)=.false. 
1667:  sp(:)=.false. 
1668:  sp2(:)=.false. 
1669:  ! build distance matrix between all C atoms and all other atoms 
1670:  do i=1,natoms 
1671:   write(labelchar(i),'(a1)') ih(m04+i-1) 
1672:  end do 
1673:  do i=1,natoms 
1674:   do j=1,natoms 
1675:    if(labelchar(i)=='C'.or.labelchar(j)=='C'.or.labelchar(i)=='N'.or.labelchar(j)=='N') then  
1676: ! caution!!! Needs further implementation if there's Calcium in the system.  
1677:     distancematrix(i,j)=sqrt((coords(3*i-2)-coords(3*j-2))**2+(coords(3*i-1)-coords(3*j-1))**2+(coords(3*i)-coords(3*j))**2) 
1678:    end if 
1679:   end do 
1680:  end do 
1681:  write(myunitnew,*) labelchar(:) 
1682:  ! loop through all atoms, check label, and find the four neighbours for the C or N atom closer than 1.7 A 
1683:  atomindex(:,:)=0 
1684:  do i=1,natoms 
1685:   if(labelchar(i)=='C'.or.labelchar(i)=='N') then 
1686:    i1=1 
1687:    do i2=1,4 
1688:     dist=1.0D10 
1689:       if(i1>4) exit 
1690:     do j=1,natoms 
1691:      if(i==j) cycle 
1692:      if(distancematrix(i,j)<1.7D0) then  ! there should be at most 4 atoms within this distance 
1693:       atomindex(i,i1)=j 
1694:       i1=i1+1 
1695:       distancematrix(i,j)=1.0D10 
1696:      end if 
1697:     end do 
1698:    end do 
1699:    ! now the C or N atoms have their four nearest neighbours determined. Detect non-sp3 carbons. 
1700: !        write(*,*) 'atomindex(:,:)', atomindex(i,:) 
1701:    if(labelchar(i)=='C') then 
1702:     do i1=1,4 
1703:      if(atomindex(i,i1)==0) then 
1704:       write(myunitnew,*) 'non-sp3 carbon detected, excluding from chirality check ', i 
1705:       sp2(i)=.true. 
1706:       atomindex(i,i1)=i 
1707:       exclude(i)=.true. 
1708:       ! see whether it is an sp2 or sp carbon 
1709: !      do i2=i1+1,4 
1710: !       if(atomindex(i,i2)==0) then 
1711: !!        write(myunitnew,*) 'sp carbon detected ', i 
1712: !        sp(i)=.true. 
1713: !        sp2(i)=.false. 
1714: !        exclude(i)=.true. 
1715: !       end if 
1716: !      end do 
1717: !      if(sp2(i)) write(myunitnew,'(A,I4)') 'sp2 carbon detected, excluding from chirality check ', i 
1718:      end if 
1719:     end do 
1720:    else ! if labelchar(i)='N' 
1721:     do i1=1,4 
1722:      if(atomindex(i,i1)==0) then 
1723:       exclude(i)=.true. 
1724:       atomindex(i,i1)=i 
1725:       do i2=i1+1,4 
1726:        if(atomindex(i,i2)==0) then 
1727:         exclude(i)=.true. 
1728: !        write(myunitnew,'(A,I4)') 'double-bonded N detected, excluding from chirality check ', i 
1729:         do i3=i2+1,4 
1730:          if(atomindex(i,i3)==0) then 
1731: !          write(myunitnew,'(A,I4)') 'triple-bonded N detected, excluding from chirality check ', i 
1732:          end if 
1733:         end do 
1734:        end if 
1735:       end do 
1736:      end if 
1737:     end do 
1738:    end if 
1739: !  now exclude C and N atoms which have two or more hydrogens attached 
1740: !  Warning: this does not take into account two methyl groups for example, or O in carboxylate groups.  
1741: !  To tackle this issue, perhaps a manual list of atom indices could be read in 
1742: !  at start with a new keyword. Those indices will be then excluded from the chirality check. 
1743:    hcount=0 
1744:    do i1=1,4 
1745:     if(labelchar(atomindex(i,i1))=='H') hcount=hcount+1 
1746:    end do 
1747:    if(hcount>1) exclude(i)=.true. 
1748:   else 
1749:    exclude(i)=.true. 
1750:   end if 
1751:  end do 
1752: end if ! if(init) 
1753: ! now calculate dihedrals for all carbons that are non sp 
1754: ! write(*,*) 'exclude ', exclude(:) 
1755: dihedral(:)=0.0D0 
1756: do i=1,natoms 
1757:  if(.not.exclude(i)) then 
1758: !  write(myunitnew,*) 'calling calc_dihedral', i, atomindex(i,:) 
1759:   call calc_dihedral(atomindex(i,:),dihedral(i),coords,natoms) 
1760:  end if 
1761: end do 
1762: ! save dihedrals calculated for the starting structure 
1763: if(init) then 
1764:  dihedralsave(:)=dihedral(:) 
1765: else ! if(.not.init) 
1766:  
1767: ! check whether the new dihedrals have opposite sign 
1768: ! (then we have an inversion of the chiral centre for sp3 carbon and cis-trans isomerisation for sp2) 
1769:  do i=1,natoms 
1770:   if(.not.exclude(i)) then 
1771:    if(dihedral(i)/dihedralsave(i)<0) then  
1772:     goodstructure=.false. 
1773:     write(myunitnew,*) 'chirality inversion (or cis-trans isomerisation) detected for atom ', i,dihedral(i),dihedralsave(i) 
1774:    end if 
1775:   end if 
1776:  end do 
1777: end if 
1778: end subroutine check_chirality_generic 
1779:  
1780: subroutine check_chirality(coords,natoms,goodstructure)1604: subroutine check_chirality(coords,natoms,goodstructure)
1781: use modamber9,only : myunitnew,i02,m02,m04,ih,ix,uachiral1605: use modamber9,only : myunitnew,i02,m02,m04,ih,ix,uachiral
1782: implicit none1606: implicit none
1783: 1607: 
1784: integer,intent(in) :: natoms1608: integer,intent(in) :: natoms
1785: double precision,intent(in) :: coords(3*natoms)1609: double precision,intent(in) :: coords(3*natoms)
1786: integer :: i,i1,j,dihed(3,4),currentresidue,resnumber1610: integer :: i,i1,j,dihed(3,4),currentresidue,resnumber
1787: double precision :: angle(4)1611: double precision :: angle(4)
1788: logical :: goodstructure1612: logical :: goodstructure
1789: 1613: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0