hdiff output

r31203/amberinterface.f 2016-09-26 17:30:13.593893669 +0100 r31202/amberinterface.f 2016-09-26 17:30:13.857897199 +0100
649: E_IGB = ener(26)649: E_IGB = ener(26)
650: E_BOND = ener(27)650: E_BOND = ener(27)
651: E_ANGLE = ener(28)651: E_ANGLE = ener(28)
652: E_DIHEDRAL = ener(29)652: E_DIHEDRAL = ener(29)
653: E_VDW = ener(24)653: E_VDW = ener(24)
654: E_14_VDW = ener(30)654: E_14_VDW = ener(30)
655: E_ELEC = ener(25)655: E_ELEC = ener(25)
656: E_14_ELEC = ener(31)656: E_14_ELEC = ener(31)
657: AMBER_IGB = igb657: AMBER_IGB = igb
658: 658: 
659: !! catch cold fusion659: ! catch cold fusion
660: !IF(.NOT.(ereal>-1.0D20 .AND. ereal<1.0D20)) THEN660: IF(.NOT.(ereal>-1.0D20 .AND. ereal<1.0D20)) THEN
661: !        write(MYUNITNEW,*) 'sf344> in amberenergies: cold fusion detected ',ereal!,natom,y(1:3*natom)661:         write(MYUNITNEW,*) 'sf344> in amberenergies: cold fusion detected ',ereal!,natom,y(1:3*natom)
662: !!     call checkdistances662: !     call checkdistances
663: !        ereal = 1.0D7663:         ereal = 1.0D7
664: !        grad(1:3*natom) = 1.0D0664:         grad(1:3*natom) = 1.0D0
665: !END IF665: END IF
666: 666: 
667: END SUBROUTINE AMBERENERGIES667: END SUBROUTINE AMBERENERGIES
668: 668: 
669: 669: 
670: SUBROUTINE TAKESTEPAMBER(JP,y,movableatomlist1,nmovableatoms,ligmovet,mdstept1,randomseedt2,blockmovet,nblocks,atomsinblock1)670: SUBROUTINE TAKESTEPAMBER(JP,y,movableatomlist,nmovableatoms,ligmovet,mdstept1,randomseedt2,blockmovet,nblocks,atomsinblock)
671: !671: !
672: ! Takestep routine for the AMBER interface.  672: ! Takestep routine for the AMBER interface.  
673: ! Short MD steps are being taken, controlled by min_md.in, followed by random xyz moves.673: ! Short MD steps are being taken, controlled by min_md.in, followed by random xyz moves.
674: ! The random xyz moves are being taken in subroutine TAKESTEP, which is being called 674: ! The random xyz moves are being taken in subroutine TAKESTEP, which is being called 
675: ! after this routine.675: ! after this routine.
676: !676: !
677: ! LIGAND ROTATION is now implemented as well. This seems to be way more efficient 677: ! LIGAND ROTATION is now implemented as well. This seems to be way more efficient 
678: !  than random xyz moves or MD runs to sample the energy landscape of a docked ligand + protein. 678: !  than random xyz moves or MD runs to sample the energy landscape of a docked ligand + protein. 
679: !  The atoms of the ligand have to be specified in the file 'movableatoms', the keywords679: !  The atoms of the ligand have to be specified in the file 'movableatoms', the keywords
680: !  'MOVABLEATOMS' and 'LIGMOVE' also need to be present in the data file.680: !  'MOVABLEATOMS' and 'LIGMOVE' also need to be present in the data file.
681: !  If one wants to explore the ligand configuration space only by rotation, it is advised681: !  If one wants to explore the ligand configuration space only by rotation, it is advised
682: !  to keep the stepsize small (STEP 0.01) and use 'AMBERMDSTEPS' as well, setting the number682: !  to keep the stepsize small (STEP 0.01) and use 'AMBERMDSTEPS' as well, setting the number
683: !  of MD steps in the min_md.in file to 1, and ntr/ntwx also to 1. In this way each starting 683: !  of MD steps in the min_md.in file to 1, and ntr/ntwx also to 1. In this way each starting 
684: !  configuration is being dumped in the coords.mdcrd file, which can be useful if one wants to 684: !  configuration is being dumped in the coords.mdcrd file, which can be useful if one wants to 
685: !  check how extensive the movement of the ligand is in the pocket.685: !  check how extensive the movement of the ligand is in the pocket.
686: use modamber9686: use modamber9
687: use commons, only : natoms, change_temp, newres_temp, perct, perccut687: use commons, only : natoms, change_temp, newres_temp
688: use nblist, only : nbflag, skinnb688: use nblist, only : nbflag, skinnb
689: use porfuncs689: use porfuncs
690: 690: 
691: implicit none691: implicit none
692: 692: 
693: 693: 
694: integer JP, itime1, now(3), nmovableatoms,i,j,movableatomlist1(nmovableatoms),nblocks,atomsinblock1(nblocks),offset1,k694: integer JP, itime1, now(3), nmovableatoms,i,j,movableatomlist(nmovableatoms),nblocks,atomsinblock(nblocks),offset1,k
695: double precision        :: y(3*natoms), grad(3*natoms), ereal,ligandcentre(3),ligandcoords(3),ligandcoordsrotated(3)695: double precision        :: y(3*natoms), grad(3*natoms), ereal,ligandcentre(3),ligandcoords(3),ligandcoordsrotated(3)
696: double precision        :: twopi,pi,DPRAND,randomphi,randompsi,randomtheta,DIST, dummyz696: double precision        :: rotationmatrix(nblocks,3,3),twopi,pi,DPRAND,randomphi,randompsi,randomtheta,DIST, dummyz
697: double precision        :: st,ct,sph,cph,sps,cps, VECBAX, VECBAY, VECBAZ, randomx, randomy, randomz, dummyx, dummyy697: double precision        :: st,ct,sph,cph,sps,cps, VECBAX, VECBAY, VECBAZ, randomx, randomy, randomz, dummyx, dummyy
698: double precision        :: distancematrix(natoms,natoms), ysave(3*natoms),cartstepsave,transstepsave698: double precision        :: atomblockcentre(nblocks,3),distancematrix(natoms,natoms)
699: logical                 :: ligmovet,mdstept1,randomseedt2,reschangedt(nres),blockmovet,overlapt699: logical                 :: ligmovet,mdstept1,randomseedt2,reschangedt(nres),blockmovet,overlapt
700: logical                 :: includedatom(natoms) 
701: character(len=10)       :: datechar,timechar,zonechar700: character(len=10)       :: datechar,timechar,zonechar
702: integer                 :: values(8),iostatus,restemp,i1,currentresidue,resnumber,overlaparray(natoms,natoms)701: integer                 :: values(8),iostatus,restemp,i1,currentresidue,resnumber,overlaparray(natoms,natoms)
703: integer                 :: overlaparraysave(natoms,natoms),nretries702: integer                 :: overlaparraysave(natoms,natoms)
704: character(len=10)       :: rotmaxchangestr,rotpselectstr,rotcutoffstr,rotcentrestr,rotoccuwstr703: character(len=10)       :: rotmaxchangestr,rotpselectstr,rotcutoffstr,rotcentrestr,rotoccuwstr
705: integer, allocatable    :: movableatomlist(:), atomsinblock(:)  
706: double precision, allocatable :: atomblockcentre(:,:), rotationmatrix(:,:,:) 
707:  
708:  
709: ! sf344> for compatibility with dynamically determining macroion nearest neighbours and moving them with LIGMOVE 
710: if(allocated(movableatomlist)) deallocate(movableatomlist) 
711: allocate(movableatomlist(nmovableatoms)) 
712: if(allocated(atomsinblock)) deallocate(atomsinblock) 
713: allocate(atomsinblock(nblocks)) 
714: if(allocated(rotationmatrix)) deallocate(rotationmatrix) 
715: allocate(rotationmatrix(nblocks,3,3)) 
716: if(allocated(atomblockcentre)) deallocate(atomblockcentre) 
717: allocate(atomblockcentre(nblocks,3)) 
718: movableatomlist(:)=movableatomlist1(:) 
719: atomsinblock(:)=atomsinblock1(:) 
720:  
721: cartstepsave=ligcartstep 
722: transstepsave=ligtransstep 
723:  
724:  
725: ! parameters704: ! parameters
726: pi=ATAN(1.0D0)*4705: pi=ATAN(1.0D0)*4
727: twopi=2.0D0*pi706: twopi=2.0D0*pi
728: ! increments with each takestepamber call707: ! increments with each takestepamber call
729: n_amb_calls=n_amb_calls+1708: n_amb_calls=n_amb_calls+1
730: 709: if(nmovableatoms==0.and.(ligmovet.or.blockmovet)) then 
731:  
732: if(.not.macroiont) then ! we determine the movable atom list dynamically for the macroion model 
733:  if(nmovableatoms==0.and.(ligmovet.or.blockmovet)) then  
734:         write(MYUNITNEW,*)&710:         write(MYUNITNEW,*)&
735:         &        ' takestepamber> ERROR - need to specify movable atoms and use the MOVABLEATOMS keyword!'711:         &        ' takestepamber> ERROR - need to specify movable atoms and use the MOVABLEATOMS keyword!'
736:         STOP712:         STOP
737:  end if713: end if
738:  if(blockmovet.and..not.ligmovet) then 714: if(blockmovet.and..not.ligmovet) then 
739:         write(MYUNITNEW,*)&715:         write(MYUNITNEW,*)&
740:         &        ' takestepamber> ERROR - BLOCKMOVE works only in conjunction with the LIGMOVE keyword!'716:         &        ' takestepamber> ERROR - BLOCKMOVE works only in conjunction with the LIGMOVE keyword!'
741:         STOP717:         STOP
742:  end if 
743: end if718: end if
744: 719: 
745: IF (STEEREDMINT) THEN720: IF (STEEREDMINT) THEN
746: 721: 
747: !722: !
748: ! (1) Move ligand out along direction from B along the vector from atom B to atom A so723: ! (1) Move ligand out along direction from B along the vector from atom B to atom A so
749: !     that the distance between SMINATOMA and SMINATOMB is SMINDISTSTART.724: !     that the distance between SMINATOMA and SMINATOMB is SMINDISTSTART.
750: ! (2) Rotate ligand (fragment A) randomly about centre of coordinates.725: ! (2) Rotate ligand (fragment A) randomly about centre of coordinates.
751: !726: !
752: ! All done below....727: ! All done below....
786:         CALL AMBERINTERFACE(natoms,4,inpcrd,MYUNITNEW)           ! retain original settings for the energy calculations761:         CALL AMBERINTERFACE(natoms,4,inpcrd,MYUNITNEW)           ! retain original settings for the energy calculations
787:         x(lcrd:lcrd+3*natom-1) = y(1:3*natom)762:         x(lcrd:lcrd+3*natom-1) = y(1:3*natom)
788: END IF763: END IF
789: 764: 
790: ! LIGAND MOVES START HERE 765: ! LIGAND MOVES START HERE 
791: ! At the moment we have rigid rotation and cartesian perturbation. There is also a translation for use in766: ! At the moment we have rigid rotation and cartesian perturbation. There is also a translation for use in
792: ! steered minimisation767: ! steered minimisation
793: 768: 
794: ! 1) RIGID BODY ROTATION769: ! 1) RIGID BODY ROTATION
795: 770: 
796: ! sf344> work out for binary systems (macroion model) the bound counterions to each macroion, and move them cooperatively 
797: !  with BLOCKMOVE. 
798:  
799:  
800:  
801: if(MACROIONT) then 
802: ! Macroions are the first in the atom list. 
803:           nblocks=nmacroions 
804:           if(allocated(movableatomlist)) deallocate(movableatomlist) 
805:           if(allocated(atomsinblock)) deallocate(atomsinblock) 
806:           if(allocated(atomblockcentre)) deallocate(atomblockcentre) 
807:           if(allocated(rotationmatrix)) deallocate(rotationmatrix) 
808:           allocate(movableatomlist(natoms),atomsinblock(nblocks),atomblockcentre(nblocks,3),rotationmatrix(nblocks,3,3)) 
809:           movableatomlist(:)=0 
810:           atomsinblock(:)=0 
811:           atomblockcentre(:,:)=0.0D0 
812:           rotationmatrix(:,:,:)=0.0D0 
813:           offset1=1 
814:           includedatom(1:natoms)=.false. 
815:           do i=1,nmacroions 
816:                movableatomlist(offset1)=i 
817:                offset1=offset1+1 
818:                atomsinblock(i)=1 
819:              do j=nmacroions+1,natoms 
820:                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) 
821:                if(distancematrix(i,j)<macroiondist.and..not.includedatom(j))then 
822:                 includedatom(j)=.true. ! include one counterion in only one block! 
823:                 movableatomlist(offset1)=j 
824:                 offset1=offset1+1 
825:                 atomsinblock(i)=atomsinblock(i)+1 
826:                end if 
827:              end do 
828:           end do 
829: !          write(MYUNITNEW,*) 'dynamically determining bound counterions to each macroion and moving them together' 
830: !          write(MYUNITNEW,*) 'natoms, nblocks', natoms, nblocks 
831: !          write(MYUNITNEW,*) 'atomsinblock(:)', atomsinblock(:) 
832: !          write(MYUNITNEW,*) 'movableatomlist(:)', movableatomlist(:) 
833: end if 
834:  
835: IF((LIGMOVET.OR.STEEREDMINT).AND.doligmove) THEN771: IF((LIGMOVET.OR.STEEREDMINT).AND.doligmove) THEN
836: 772: 
837:         IF (STEEREDMINT) WRITE(MYUNITNEW,'(A)') 'takestepamber> rotating and translating the ligand'773:         IF (STEEREDMINT) WRITE(MYUNITNEW,'(A)') 'takestepamber> rotating and translating the ligand'
838:   overlapt=.true.774:   overlapt=.true.
839:   overlaparraysave(:,:)=0775:   overlaparraysave(:,:)=0
840:   overlaparray(:,:)=0776:   overlaparray(:,:)=0
841:   distancematrix(:,:)=0.0D0777:   distancematrix(:,:)=0.0D0
842:   ysave(:)=y(:) 
843:         IF(BLOCKMOVET) THEN778:         IF(BLOCKMOVET) THEN
844:          ! Try to avoid moves that cause extremely high initial forces (or cold fusion),779:          ! Try to avoid moves that cause extremely high initial forces (or cold fusion),
845:          ! by redoing the rotation moves whenever two atoms get closer than 0.5 Angstroms.780:          ! by redoing the rotation moves whenever two atoms get closer than 0.5 Angstroms.
846:           do i=1,natoms781:           do i=1,natoms
847:            do j=1,natoms782:            do j=1,natoms
848:               if(i==j) cycle783:               if(i==j) cycle
849:               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)784:               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)
850:               if((MACROIONT.and.distancematrix(i,j)<1.0D0).or.distancematrix(i,j)<0.7D0) then 785:               if(distancematrix(i,j)<0.7D0) then 
851:                 overlaparraysave(i,j)=1  !there's an initial overlap786:                 overlaparraysave(i,j)=1  !there's an initial overlap
852: !                write(*,*) 'atoms close together before rotating ', i,j,distancematrix(i,j)787: !                write(*,*) 'atoms close together before rotating ', i,j,distancematrix(i,j)
853:               end if788:               end if
854:            end do789:            end do
855:           end do790:           end do
856:         END IF791:         END IF
857:  
858:  nretries=0 
859:  DO WHILE(OVERLAPT)792:  DO WHILE(OVERLAPT)
860:    IF(nretries>10) then ! scale ligcartstep and ligtransstep 
861:       nretries=0 
862:       ligcartstep=ligcartstep*0.9 
863:       ligtransstep=ligtransstep*0.9 
864: !      WRITE(MYUNITNEW,*) 'scaling down ligcartstep and ligtransstep', ligcartstep, ligtransstep 
865:    END IF 
866:      overlapt=.false.793:      overlapt=.false.
867:      DO i=1,nblocks ! nblocks=1 if BLOCKMOVET=.FALSE. 794:      DO i=1,nblocks ! nblocks=1 if BLOCKMOVET=.FALSE. 
868:         randomphi=(DPRAND()-0.5)*twopi*ligrotscale795:         randomphi=(DPRAND()-0.5)*twopi*ligrotscale
869:         randomtheta=(DPRAND()-0.5)*pi*ligrotscale796:         randomtheta=(DPRAND()-0.5)*pi*ligrotscale
870:         randompsi=(DPRAND()-0.5)*twopi*ligrotscale797:         randompsi=(DPRAND()-0.5)*twopi*ligrotscale
871:         st=sin(randomtheta)798:         st=sin(randomtheta)
872:         ct=cos(randomtheta)799:         ct=cos(randomtheta)
873:         sph=sin(randomphi)800:         sph=sin(randomphi)
874:         cph=cos(randomphi)801:         cph=cos(randomphi)
875:         sps=sin(randompsi)802:         sps=sin(randompsi)
885:         rotationmatrix(i,1,3)=st*sph812:         rotationmatrix(i,1,3)=st*sph
886:         rotationmatrix(i,2,3)=-st*cph813:         rotationmatrix(i,2,3)=-st*cph
887:         rotationmatrix(i,3,3)=ct814:         rotationmatrix(i,3,3)=ct
888: !        WRITE(*,*) 'rotation matrix: ', rotationmatrix(i,:,:)815: !        WRITE(*,*) 'rotation matrix: ', rotationmatrix(i,:,:)
889:      END DO816:      END DO
890: !         y(1:3*natom)=coords1(1:3*natom)817: !         y(1:3*natom)=coords1(1:3*natom)
891: 818: 
892: !        y(1:3*natom)=x(lcrd:lcrd+3*natom-1)819: !        y(1:3*natom)=x(lcrd:lcrd+3*natom-1)
893: 820: 
894: ! work out the centre of coordinates for the ligand821: ! work out the centre of coordinates for the ligand
895:         IF(BLOCKMOVET) THEN 
896:  
897: 822: 
 823:         IF(BLOCKMOVET) THEN
898:           atomblockcentre(:,:)=0.0D0824:           atomblockcentre(:,:)=0.0D0
899:           offset1=0825:           offset1=0
900: !          overlapt=.false.826: !          overlapt=.false.
901:           do i=1,nblocks827:           do i=1,nblocks
902:              if(i>1) offset1=offset1+atomsinblock(i-1)828:              if(i>1) offset1=offset1+atomsinblock(i-1)
903:              do k=1,atomsinblock(i)829:              do k=1,atomsinblock(i)
904:                 j=movableatomlist(k+offset1)830:                 j=movableatomlist(k+offset1)
905: !                WRITE(*,*) k, j831: !                WRITE(*,*) k, j
906:                 atomblockcentre(i,1)=atomblockcentre(i,1)+y(3*j-2)832:                 atomblockcentre(i,1)=atomblockcentre(i,1)+y(3*j-2)
907:                 atomblockcentre(i,2)=atomblockcentre(i,2)+y(3*j-1)833:                 atomblockcentre(i,2)=atomblockcentre(i,2)+y(3*j-1)
933:               else859:               else
934:                  overlaparray(i,j)=0 860:                  overlaparray(i,j)=0 
935:               end if861:               end if
936:               if(overlaparraysave(i,j)-overlaparray(i,j)<0) then 862:               if(overlaparraysave(i,j)-overlaparray(i,j)<0) then 
937:                 overlapt=.true.863:                 overlapt=.true.
938: !                write(*,*) 'atoms close together after rotating ', i,j,distancematrix(i,j)864: !                write(*,*) 'atoms close together after rotating ', i,j,distancematrix(i,j)
939:                 GOTO 95865:                 GOTO 95
940:               end if866:               end if
941:            end do867:            end do
942:          end do868:          end do
943:  
944:         ELSE869:         ELSE
945:           ligandcentre(:)=0.0D0870:           ligandcentre(:)=0.0D0
946:           do i=1,nmovableatoms871:           do i=1,nmovableatoms
947:                 j=movableatomlist(i)872:                 j=movableatomlist(i)
948:                 ligandcentre(1)=ligandcentre(1)+y(3*j-2)873:                 ligandcentre(1)=ligandcentre(1)+y(3*j-2)
949:                 ligandcentre(2)=ligandcentre(2)+y(3*j-1)874:                 ligandcentre(2)=ligandcentre(2)+y(3*j-1)
950:                 ligandcentre(3)=ligandcentre(3)+y(3*j  )875:                 ligandcentre(3)=ligandcentre(3)+y(3*j  )
951:           end do876:           end do
952:                 ligandcentre(:) = ligandcentre(:)/nmovableatoms877:                 ligandcentre(:) = ligandcentre(:)/nmovableatoms
953: !        write(MYUNITNEW,*) 'ligandcentre = ', ligandcentre(:)878: !        write(MYUNITNEW,*) 'ligandcentre = ', ligandcentre(:)
961:                 ligandcoordsrotated=MATMUL(rotationmatrix(1,:,:),ligandcoords)886:                 ligandcoordsrotated=MATMUL(rotationmatrix(1,:,:),ligandcoords)
962: !          write(MYUNITNEW,*) 'rotated coords: ', ligandcoordsrotated(:)887: !          write(MYUNITNEW,*) 'rotated coords: ', ligandcoordsrotated(:)
963: !          write(MYUNITNEW,*) 'original coords: ', ligandcoords(:)888: !          write(MYUNITNEW,*) 'original coords: ', ligandcoords(:)
964: 889: 
965: ! translate the ligand atoms back to the centre of the original ligand                890: ! translate the ligand atoms back to the centre of the original ligand                
966:                 y(3*j-2)=ligandcoordsrotated(1)+ligandcentre(1)891:                 y(3*j-2)=ligandcoordsrotated(1)+ligandcentre(1)
967:                 y(3*j-1)=ligandcoordsrotated(2)+ligandcentre(2)892:                 y(3*j-1)=ligandcoordsrotated(2)+ligandcentre(2)
968:                 y(3*j  )=ligandcoordsrotated(3)+ligandcentre(3)893:                 y(3*j  )=ligandcoordsrotated(3)+ligandcentre(3)
969:           end do894:           end do
970:         END IF895:         END IF
 896: 95        CONTINUE
 897:  END DO  ! while(overlapt)
971: ! 2) RANDOM CARTESIAN PERTURBATION898: ! 2) RANDOM CARTESIAN PERTURBATION
972:         if ((LIGMOVET).and.(ligcartstep.gt.0)) then899:         if ((LIGMOVET).and.(ligcartstep.gt.0)) then
973:            do i=1,nmovableatoms900:            do i=1,nmovableatoms
974:               j=movableatomlist(i)901:               j=movableatomlist(i)
975:               randomx=DPRAND()902:               randomx=DPRAND()
976:               randomy=DPRAND()903:               randomy=DPRAND()
977:               randomz=DPRAND()904:               randomz=DPRAND()
978:               y(3*j-2)=y(3*j-2)+2*randomx*ligcartstep-ligcartstep905:               y(3*j-2)=y(3*j-2)+2*randomx*ligcartstep-ligcartstep
979:               y(3*j-1)=y(3*j-1)+2*randomy*ligcartstep-ligcartstep906:               y(3*j-1)=y(3*j-1)+2*randomy*ligcartstep-ligcartstep
980:               y(3*j  )=y(3*j  )+2*randomz*ligcartstep-ligcartstep907:               y(3*j  )=y(3*j  )+2*randomz*ligcartstep-ligcartstep
987:            randomz=DPRAND()914:            randomz=DPRAND()
988:            do i=1,nmovableatoms915:            do i=1,nmovableatoms
989:               j=movableatomlist(i)916:               j=movableatomlist(i)
990:               y(3*j-2)=y(3*j-2)+2*randomx*ligtransstep-ligtransstep917:               y(3*j-2)=y(3*j-2)+2*randomx*ligtransstep-ligtransstep
991:               y(3*j-1)=y(3*j-1)+2*randomy*ligtransstep-ligtransstep918:               y(3*j-1)=y(3*j-1)+2*randomy*ligtransstep-ligtransstep
992:               y(3*j  )=y(3*j  )+2*randomz*ligtransstep-ligtransstep919:               y(3*j  )=y(3*j  )+2*randomz*ligtransstep-ligtransstep
993:            enddo920:            enddo
994:         endif921:         endif
995:         if ((LIGMOVET).and.(ligtransstep.gt.0).and.BLOCKMOVET) then922:         if ((LIGMOVET).and.(ligtransstep.gt.0).and.BLOCKMOVET) then
996:           offset1=0923:           offset1=0
997: !          WRITE(*,*) 'nblocks ', nblocks 
998:           do i=1,nblocks924:           do i=1,nblocks
999:              if(i>1) offset1=offset1+atomsinblock(i-1)925:              if(i>1) offset1=offset1+atomsinblock(i-1)
1000:              randomx=DPRAND()926:              randomx=DPRAND()
1001:              randomy=DPRAND()927:              randomy=DPRAND()
1002:              randomz=DPRAND()928:              randomz=DPRAND()
1003:              do k=1,atomsinblock(i)929:              do k=1,atomsinblock(i)
1004:                 j=movableatomlist(k+offset1)930:                 j=movableatomlist(k+offset1)
1005:                 y(3*j-2)=y(3*j-2)+2*randomx*ligtransstep-ligtransstep931:                 y(3*j-2)=y(3*j-2)+2*randomx*ligtransstep-ligtransstep
1006:                 y(3*j-1)=y(3*j-1)+2*randomy*ligtransstep-ligtransstep932:                 y(3*j-1)=y(3*j-1)+2*randomy*ligtransstep-ligtransstep
1007:                 y(3*j  )=y(3*j  )+2*randomz*ligtransstep-ligtransstep933:                 y(3*j  )=y(3*j  )+2*randomz*ligtransstep-ligtransstep
1008:              end do934:              end do
1009:           end do935:           end do
1010:         end if936:         end if
1011:  ! check for percolation, and take a new step until we start off from a percolated structure937:   
1012:         IF(MACROIONT) THEN 
1013:               PERCT=.TRUE. 
1014:               CALL PERC(y,NATOMS,PERCCUT,PERCT,.FALSE.,MYUNITNEW,.FALSE.) 
1015: !              WRITE(*,*) 'perccut=', perccut 
1016:               IF(.NOT.PERCT) then 
1017: !               write(MYUNITNEW,*) ' takestepamber> disconnected structure detected, undoing step' 
1018:                 y(:)=ysave(:) 
1019:                 overlapt=.true. 
1020:                 nretries=nretries+1  ! make sure we are not getting stuck in an infinite loop 
1021:                 GOTO 95 
1022:               ELSE 
1023:                 overlapt=.false. 
1024:               END IF 
1025:         END IF 
1026:   
1027: 95        CONTINUE 
1028:  END DO  ! while(overlapt) 
1029: ! reset ligcartstep and ligtransstep 
1030: ligcartstep=cartstepsave 
1031: ligtransstep=transstepsave 
1032: 938: 
1033: ! 4) TRANSLATION MOVE FOR STEERED MINIMISATION (rotation done above)939: ! 4) TRANSLATION MOVE FOR STEERED MINIMISATION (rotation done above)
1034: 940: 
1035: 941: 
1036: ! some sanity checks!942: ! some sanity checks!
1037: ! check atom A is in the ligand (specified in movableatoms)943: ! check atom A is in the ligand (specified in movableatoms)
1038:    IF (STEEREDMINT) THEN944:    IF (STEEREDMINT) THEN
1039:       do i=1,nmovableatoms945:       do i=1,nmovableatoms
1040:          j=movableatomlist(i)946:          j=movableatomlist(i)
1041:          IF (J.EQ.SMINATOMA) GOTO 678947:          IF (J.EQ.SMINATOMA) GOTO 678
1228: !IF(nunit/=20 .and. nunit/=25) THEN1134: !IF(nunit/=20 .and. nunit/=25) THEN
1229: !        write(MYUNITNEW,*) 'AMBERfinalio> nunit has to be 20 (OPTIM) or 25 (GMIN), &1135: !        write(MYUNITNEW,*) 'AMBERfinalio> nunit has to be 20 (OPTIM) or 25 (GMIN), &
1230: !&                        using nunit=20'1136: !&                        using nunit=20'
1231: !        nunit=201137: !        nunit=20
1232: !END IF1138: !END IF
1233: 1139: 
1234: restrt1 = restrt1140: restrt1 = restrt
1235: !write(MYUNITNEW,*) 'NATOM', natom1141: !write(MYUNITNEW,*) 'NATOM', natom
1236: 1142: 
1237:     x(lcrd:lcrd+3*natom-1) = coords2(1:3*natom)1143:     x(lcrd:lcrd+3*natom-1) = coords2(1:3*natom)
1238: IF(FILTH==0) THEN      ! dump output files for GMIN1144: IF(nunit==25) THEN      ! dump output files for GMIN
1239:     write(pdbfile,'(A)') 'lowest'//trim(adjustl(nsavestring))//'.pdb'1145:     write(pdbfile,'(A)') 'lowest'//trim(adjustl(nsavestring))//'.pdb'
1240:     write(restrt2,'(A)') 'min'//trim(adjustl(nsavestring))//'.rst'1146:     write(restrt2,'(A)') 'min'//trim(adjustl(nsavestring))//'.rst'
1241: ELSE IF(FILTH==1) THEN ! GMIN with MPI1147: ELSE IF(nunit/=20) THEN ! GMIN with MPI
1242:     write(nodestring,'(I6)') node ! MYNODE in main.F1148:     write(nodestring,'(I6)') node ! MYNODE in main.F
1243:     write(pdbfile,'(A)') 'lowest'//trim(adjustl(nsavestring))//'.'//trim(adjustl(nodestring))//'.pdb'1149:     write(pdbfile,'(A)') 'lowest'//trim(adjustl(nsavestring))//'.'//trim(adjustl(nodestring))//'.pdb'
1244:     write(restrt2,'(A)') 'min'//trim(adjustl(nsavestring))//'.'//trim(adjustl(nodestring))//'.rst'1150:     write(restrt2,'(A)') 'min'//trim(adjustl(nsavestring))//'.'//trim(adjustl(nodestring))//'.rst'
1245: ELSE IF(FILTH==2) THEN1151: ELSE                    ! nunit=20, dump files for OPTIM
 1152:  IF(FILTH.EQ.0) THEN
1246:     open(unit=nunit,file='points.final.xyz',status='unknown')1153:     open(unit=nunit,file='points.final.xyz',status='unknown')
1247:     write(nunit,'(I6)') natom1154:     write(nunit,'(I6)') natom
1248:     write(nunit,'(A)') 'points.final'1155:     write(nunit,'(A)') 'points.final'
1249:     write(pdbfile,'(A)') 'points.final.pdb'1156:     write(pdbfile,'(A)') 'points.final.pdb'
1250:     write(restrt2,'(A)') 'points.final.rst'1157:     write(restrt2,'(A)') 'points.final.rst'
1251: ELSE1158:  ELSE
1252:     open(unit=nunit,file='points.final.xyz.'//trim(adjustl(ostring)),status='unknown')1159:     open(unit=nunit,file='points.final.xyz.'//trim(adjustl(ostring)),status='unknown')
1253:     write(nunit,'(I6)') natom1160:     write(nunit,'(I6)') natom
1254:     write(nunit,'(A)') 'points.final.'//trim(adjustl(ostring))1161:     write(nunit,'(A)') 'points.final.'//trim(adjustl(ostring))
1255:     write(pdbfile,'(A)') 'points.final.pdb.'//trim(adjustl(ostring))1162:     write(pdbfile,'(A)') 'points.final.pdb.'//trim(adjustl(ostring))
1256:     write(restrt2,'(A)') 'points.final.rst.'//trim(adjustl(ostring))1163:     write(restrt2,'(A)') 'points.final.rst.'//trim(adjustl(ostring))
 1164:  END IF
1257: END IF1165: END IF
1258: 1166: 
1259:     restrt = restrt21167:     restrt = restrt2
1260: !WRITE(nunit,*) 'nsave=', nsave, trim(adjustl(pdbfile))1168: !WRITE(nunit,*) 'nsave=', nsave, trim(adjustl(pdbfile))
1261: !WRITE(nunit,*) coords2(1:3)1169: !WRITE(nunit,*) coords2(1:3)
1262: !WRITE(nunit,*) x(lcrd:lcrd+2)1170: !WRITE(nunit,*) x(lcrd:lcrd+2)
1263: DO i1=1,natom1171: DO i1=1,natom
1264:         WRITE(nunit,'(A4,3F20.10)') ih(m04+i1-1),x(lcrd+3*i1-3),x(lcrd+3*i1-2),x(lcrd+3*i1-1) 1172:         WRITE(nunit,'(A4,3F20.10)') ih(m04+i1-1),x(lcrd+3*i1-3),x(lcrd+3*i1-2),x(lcrd+3*i1-1) 
1265: END DO1173: END DO
1266: ambpdb_unit=GETUNIT()1174: ambpdb_unit=GETUNIT()


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0