hdiff output

r22793/amberinterface.f 2017-03-30 13:30:06.104194724 +0100 r22792/amberinterface.f 2017-03-30 13:30:06.404198698 +0100
639: ! catch cold fusion:639: ! catch cold fusion:
640: IF(ereal<-1.0D10) THEN640: IF(ereal<-1.0D10) THEN
641:         write(MYUNITNEW,*) 'sf344> in amberenergies: cold fusion detected, checking distances',ereal!,natom,y(1:3*natom)641:         write(MYUNITNEW,*) 'sf344> in amberenergies: cold fusion detected, checking distances',ereal!,natom,y(1:3*natom)
642: !     call checkdistances642: !     call checkdistances
643:         ereal = -1.0D7643:         ereal = -1.0D7
644:         grad(1:3*natom) = 1.0D0644:         grad(1:3*natom) = 1.0D0
645: END IF645: END IF
646: END SUBROUTINE AMBERENERGIES646: END SUBROUTINE AMBERENERGIES
647: 647: 
648: 648: 
649: SUBROUTINE TAKESTEPAMBER(JP,y,movableatomlist,nmovableatoms,ligmovet,mdstept1,randomseedt2,blockmovet,nblocks,atomsinblock)649: SUBROUTINE TAKESTEPAMBER(JP,y,movableatomlist,nmovableatoms,ligmovet,mdstept1,randomseedt2)
650: !650: !
651: ! Takestep routine for the AMBER interface.  651: ! Takestep routine for the AMBER interface.  
652: ! Short MD steps are being taken, controlled by min_md.in, followed by random xyz moves.652: ! Short MD steps are being taken, controlled by min_md.in, followed by random xyz moves.
653: ! The random xyz moves are being taken in subroutine TAKESTEP, which is being called 653: ! The random xyz moves are being taken in subroutine TAKESTEP, which is being called 
654: ! after this routine.654: ! after this routine.
655: !655: !
656: ! LIGAND ROTATION is now implemented as well. This seems to be way more efficient 656: ! LIGAND ROTATION is now implemented as well. This seems to be way more efficient 
657: !  than random xyz moves or MD runs to sample the energy landscape of a docked ligand + protein. 657: !  than random xyz moves or MD runs to sample the energy landscape of a docked ligand + protein. 
658: !  The atoms of the ligand have to be specified in the file 'movableatoms', the keywords658: !  The atoms of the ligand have to be specified in the file 'movableatoms', the keywords
659: !  'MOVABLEATOMS' and 'LIGMOVE' also need to be present in the data file.659: !  'MOVABLEATOMS' and 'LIGMOVE' also need to be present in the data file.
663: !  configuration is being dumped in the coords.mdcrd file, which can be useful if one wants to 663: !  configuration is being dumped in the coords.mdcrd file, which can be useful if one wants to 
664: !  check how extensive the movement of the ligand is in the pocket.664: !  check how extensive the movement of the ligand is in the pocket.
665: use modamber9665: use modamber9
666: use commons, only : natoms, change_temp, newres_temp666: use commons, only : natoms, change_temp, newres_temp
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)
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(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)677: logical                 :: ligmovet,mdstept1,randomseedt2,reschangedt(nres)
678: logical                 :: ligmovet,mdstept1,randomseedt2,reschangedt(nres),blockmovet 
679: character(len=10)       :: datechar,timechar,zonechar678: character(len=10)       :: datechar,timechar,zonechar
680: integer                 :: values(8),iostatus,restemp,i1,currentresidue,resnumber679: integer                 :: values(8),iostatus,restemp,i1,currentresidue,resnumber
681: character(len=10)       :: rotmaxchangestr,rotpselectstr,rotcutoffstr,rotcentrestr,rotoccuwstr680: character(len=10)       :: rotmaxchangestr,rotpselectstr,rotcutoffstr,rotcentrestr,rotoccuwstr
682: ! parameters681: ! parameters
683: pi=ATAN(1.0D0)*4682: pi=ATAN(1.0D0)*4
684: twopi=2.0D0*pi683: twopi=2.0D0*pi
685: ! increments with each takestepamber call684: ! increments with each takestepamber call
686: n_amb_calls=n_amb_calls+1685: n_amb_calls=n_amb_calls+1
687: if(nmovableatoms==0.and.(ligmovet.or.blockmovet)) then 686: if(nmovableatoms==0.and.ligmovet) then 
688:         write(MYUNITNEW,*)&687:         write(MYUNITNEW,*)&
689:         &        ' takestepamber> ERROR - need to specify movable atoms and use the MOVABLEATOMS keyword!'688:         &        ' takestepamber> ERROR - need to specify movable atoms and use the MOVABLEATOMS keyword!'
690:         STOP689:         STOP
691: end if690: end if
692: if(blockmovet.and..not.ligmovet) then  
693:         write(MYUNITNEW,*)& 
694:         &        ' takestepamber> ERROR - BLOCKMOVE works only in conjunction with the LIGMOVE keyword!' 
695:         STOP 
696: end if 
697: 691: 
698: IF (STEEREDMINT) THEN692: IF (STEEREDMINT) THEN
699: 693: 
700: !694: !
701: ! (1) Move ligand out along direction from B along the vector from atom B to atom A so695: ! (1) Move ligand out along direction from B along the vector from atom B to atom A so
702: !     that the distance between SMINATOMA and SMINATOMB is SMINDISTSTART.696: !     that the distance between SMINATOMA and SMINATOMB is SMINDISTSTART.
703: ! (2) Rotate ligand (fragment A) randomly about centre of coordinates.697: ! (2) Rotate ligand (fragment A) randomly about centre of coordinates.
704: !698: !
705: ! All done below....699: ! All done below....
706: !700: !
742: 736: 
743: ! LIGAND MOVES START HERE 737: ! LIGAND MOVES START HERE 
744: ! At the moment we have rigid rotation and cartesian perturbation. There is also a translation for use in738: ! At the moment we have rigid rotation and cartesian perturbation. There is also a translation for use in
745: ! steered minimisation739: ! steered minimisation
746: 740: 
747: ! 1) RIGID BODY ROTATION741: ! 1) RIGID BODY ROTATION
748: 742: 
749: IF((LIGMOVET.OR.STEEREDMINT).AND.doligmove) THEN743: IF((LIGMOVET.OR.STEEREDMINT).AND.doligmove) THEN
750: 744: 
751:         IF (STEEREDMINT) WRITE(MYUNITNEW,'(A)') 'takestepamber> rotating and translating the ligand'745:         IF (STEEREDMINT) WRITE(MYUNITNEW,'(A)') 'takestepamber> rotating and translating the ligand'
752:      DO i=1,nblocks ! nblocks=1 if BLOCKMOVET=.FALSE. 746: 
753:         randomphi=(DPRAND()-0.5)*twopi*ligrotscale747:         randomphi=(DPRAND()-0.5)*twopi*ligrotscale
754:         randomtheta=(DPRAND()-0.5)*pi*ligrotscale748:         randomtheta=(DPRAND()-0.5)*pi*ligrotscale
755:         randompsi=(DPRAND()-0.5)*twopi*ligrotscale749:         randompsi=(DPRAND()-0.5)*twopi*ligrotscale
756:         st=sin(randomtheta)750:         st=sin(randomtheta)
757:         ct=cos(randomtheta)751:         ct=cos(randomtheta)
758:         sph=sin(randomphi)752:         sph=sin(randomphi)
759:         cph=cos(randomphi)753:         cph=cos(randomphi)
760:         sps=sin(randompsi)754:         sps=sin(randompsi)
761:         cps=cos(randompsi)755:         cps=cos(randompsi)
762: 756: 
763: 757: 
764:         rotationmatrix(i,1,1)=cps*cph-ct*sph*sps758:         rotationmatrix(1,1)=cps*cph-ct*sph*sps
765:         rotationmatrix(i,2,1)=cps*sph+ct*cph*sps759:         rotationmatrix(2,1)=cps*sph+ct*cph*sps
766:         rotationmatrix(i,3,1)=sps*st760:         rotationmatrix(3,1)=sps*st
767:         rotationmatrix(i,1,2)=-sps*cph-ct*sph*cps761:         rotationmatrix(1,2)=-sps*cph-ct*sph*cps
768:         rotationmatrix(i,2,2)=-sps*sph+ct*cph*cps762:         rotationmatrix(2,2)=-sps*sph+ct*cph*cps
769:         rotationmatrix(i,3,2)=cps*st763:         rotationmatrix(3,2)=cps*st
770:         rotationmatrix(i,1,3)=st*sph764:         rotationmatrix(1,3)=st*sph
771:         rotationmatrix(i,2,3)=-st*cph765:         rotationmatrix(2,3)=-st*cph
772:         rotationmatrix(i,3,3)=ct766:         rotationmatrix(3,3)=ct
773: !        WRITE(*,*) 'rotation matrix: ', rotationmatrix(i,:,:)767: 
774:      END DO 
775: !         y(1:3*natom)=coords1(1:3*natom)768: !         y(1:3*natom)=coords1(1:3*natom)
776: 769: 
777: !        y(1:3*natom)=x(lcrd:lcrd+3*natom-1)770: !        y(1:3*natom)=x(lcrd:lcrd+3*natom-1)
778: 771: 
779: ! work out the centre of coordinates for the ligand772: ! work out the centre of coordinates for the ligand
780: 773:         ligandcentre(:)=0.0D0
781:         IF(BLOCKMOVET) THEN774:         do i=1,nmovableatoms
782:           atomblockcentre(:,:)=0.0D0 
783:           offset1=0 
784:           do i=1,nblocks 
785:              if(i>1) offset1=offset1+atomsinblock(i-1) 
786:              do k=1,atomsinblock(i) 
787:                 j=movableatomlist(k+offset1) 
788: !                WRITE(*,*) k, j 
789:                 atomblockcentre(i,1)=atomblockcentre(i,1)+y(3*j-2) 
790:                 atomblockcentre(i,2)=atomblockcentre(i,2)+y(3*j-1) 
791:                 atomblockcentre(i,3)=atomblockcentre(i,3)+y(3*j  ) 
792:              end do 
793:                 atomblockcentre(i,:)=atomblockcentre(i,:)/atomsinblock(i) 
794:              do k=1,atomsinblock(i) 
795:                 j=movableatomlist(k+offset1) 
796: ! move each block to the origin 
797:                 ligandcoords(1)=y(3*j-2)-atomblockcentre(i,1) 
798:                 ligandcoords(2)=y(3*j-1)-atomblockcentre(i,2) 
799:                 ligandcoords(3)=y(3*j  )-atomblockcentre(i,3) 
800: ! rotate the block of atoms with random rotation matrix 
801:                 ligandcoordsrotated=MATMUL(rotationmatrix(i,:,:),ligandcoords) 
802: ! translate back the block of atoms to the original centre of their coordinates 
803:                 y(3*j-2)=ligandcoordsrotated(1)+atomblockcentre(i,1) 
804:                 y(3*j-1)=ligandcoordsrotated(2)+atomblockcentre(i,2) 
805:                 y(3*j  )=ligandcoordsrotated(3)+atomblockcentre(i,3) 
806:              end do             
807:           end do 
808:         ELSE 
809:           ligandcentre(:)=0.0D0 
810:           do i=1,nmovableatoms 
811:                 j=movableatomlist(i)775:                 j=movableatomlist(i)
812:                 ligandcentre(1)=ligandcentre(1)+y(3*j-2)776:                 ligandcentre(1)=ligandcentre(1)+y(3*j-2)
813:                 ligandcentre(2)=ligandcentre(2)+y(3*j-1)777:                 ligandcentre(2)=ligandcentre(2)+y(3*j-1)
814:                 ligandcentre(3)=ligandcentre(3)+y(3*j  )778:                 ligandcentre(3)=ligandcentre(3)+y(3*j  )
815:           end do779:         end do
816:                 ligandcentre(:) = ligandcentre(:)/nmovableatoms780:                 ligandcentre(:) = ligandcentre(:)/nmovableatoms
817: !        write(MYUNITNEW,*) 'ligandcentre = ', ligandcentre(:)781: !        write(MYUNITNEW,*) 'ligandcentre = ', ligandcentre(:)
818:           do i=1,nmovableatoms782:         do i=1,nmovableatoms
819: ! move the ligand atoms to the origin783: ! move the ligand atoms to the origin
820:                 j=movableatomlist(i)784:                 j=movableatomlist(i)
821:                 ligandcoords(1)=y(3*j-2)-ligandcentre(1)785:                 ligandcoords(1)=y(3*j-2)-ligandcentre(1)
822:                 ligandcoords(2)=y(3*j-1)-ligandcentre(2)786:                 ligandcoords(2)=y(3*j-1)-ligandcentre(2)
823:                 ligandcoords(3)=y(3*j  )-ligandcentre(3)787:                 ligandcoords(3)=y(3*j  )-ligandcentre(3)
824: ! rotate the ligand with the random rotation matrix788: ! rotate the ligand with the random rotation matrix
825:                 ligandcoordsrotated=MATMUL(rotationmatrix(1,:,:),ligandcoords)789:                 ligandcoordsrotated=MATMUL(rotationmatrix,ligandcoords)
826: !          write(MYUNITNEW,*) 'rotated coords: ', ligandcoordsrotated(:)790: !          write(MYUNITNEW,*) 'rotated coords: ', ligandcoordsrotated(:)
827: !          write(MYUNITNEW,*) 'original coords: ', ligandcoords(:)791: !          write(MYUNITNEW,*) 'original coords: ', ligandcoords(:)
828: 792: 
829: ! translate the ligand atoms back to the centre of the original ligand                793: ! translate the ligand atoms back to the centre of the original ligand                
830:                 y(3*j-2)=ligandcoordsrotated(1)+ligandcentre(1)794:                 y(3*j-2)=ligandcoordsrotated(1)+ligandcentre(1)
831:                 y(3*j-1)=ligandcoordsrotated(2)+ligandcentre(2)795:                 y(3*j-1)=ligandcoordsrotated(2)+ligandcentre(2)
832:                 y(3*j  )=ligandcoordsrotated(3)+ligandcentre(3)796:                 y(3*j  )=ligandcoordsrotated(3)+ligandcentre(3)
833:           end do797:         end do
834:         END IF 
835: 798: 
836: ! 2) RANDOM CARTESIAN PERTURBATION799: ! 2) RANDOM CARTESIAN PERTURBATION
837:         if ((LIGMOVET).and.(ligcartstep.gt.0)) then800:         if ((LIGMOVET).and.(ligcartstep.gt.0)) then
838:            do i=1,nmovableatoms801:            do i=1,nmovableatoms
839:               j=movableatomlist(i)802:               j=movableatomlist(i)
840:               randomx=DPRAND()803:               randomx=DPRAND()
841:               randomy=DPRAND()804:               randomy=DPRAND()
842:               randomz=DPRAND()805:               randomz=DPRAND()
843:               y(3*j-2)=y(3*j-2)+2*randomx*ligcartstep-ligcartstep806:               y(3*j-2)=y(3*j-2)+2*randomx*ligcartstep-ligcartstep
844:               y(3*j-1)=y(3*j-1)+2*randomy*ligcartstep-ligcartstep807:               y(3*j-1)=y(3*j-1)+2*randomy*ligcartstep-ligcartstep
845:               y(3*j  )=y(3*j  )+2*randomz*ligcartstep-ligcartstep808:               y(3*j  )=y(3*j  )+2*randomz*ligcartstep-ligcartstep
846:            enddo809:            enddo
847:         endif810:         endif
 811: 
848: ! 3) RIGID LIGAND TRANSLATION812: ! 3) RIGID LIGAND TRANSLATION
849:         if ((LIGMOVET).and.(ligtransstep.gt.0).and..not.BLOCKMOVET) then813:         if ((LIGMOVET).and.(ligtransstep.gt.0)) then
850:            randomx=DPRAND()814:            randomx=DPRAND()
851:            randomy=DPRAND()815:            randomy=DPRAND()
852:            randomz=DPRAND()816:            randomz=DPRAND()
853:            do i=1,nmovableatoms817:            do i=1,nmovableatoms
854:               j=movableatomlist(i)818:               j=movableatomlist(i)
855:               y(3*j-2)=y(3*j-2)+2*randomx*ligtransstep-ligtransstep819:               y(3*j-2)=y(3*j-2)+2*randomx*ligtransstep-ligtransstep
856:               y(3*j-1)=y(3*j-1)+2*randomy*ligtransstep-ligtransstep820:               y(3*j-1)=y(3*j-1)+2*randomy*ligtransstep-ligtransstep
857:               y(3*j  )=y(3*j  )+2*randomz*ligtransstep-ligtransstep821:               y(3*j  )=y(3*j  )+2*randomz*ligtransstep-ligtransstep
858:            enddo822:            enddo
859:         endif823:         endif
860:         if ((LIGMOVET).and.(ligtransstep.gt.0).and.BLOCKMOVET) then824:    
861:           offset1=0 
862:           do i=1,nblocks 
863:              if(i>1) offset1=offset1+atomsinblock(i-1) 
864:              randomx=DPRAND() 
865:              randomy=DPRAND() 
866:              randomz=DPRAND() 
867:              do k=1,atomsinblock(i) 
868:                 j=movableatomlist(k+offset1) 
869:                 y(3*j-2)=y(3*j-2)+2*randomx*ligtransstep-ligtransstep 
870:                 y(3*j-1)=y(3*j-1)+2*randomy*ligtransstep-ligtransstep 
871:                 y(3*j  )=y(3*j  )+2*randomz*ligtransstep-ligtransstep 
872:              end do 
873:           end do 
874:         end if 
875:    
876: 825: 
877: ! 4) TRANSLATION MOVE FOR STEERED MINIMISATION (rotation done above)826: ! 4) TRANSLATION MOVE FOR STEERED MINIMISATION (rotation done above)
878: 827: 
879: 828: 
880: ! some sanity checks!829: ! some sanity checks!
881: ! check atom A is in the ligand (specified in movableatoms)830: ! check atom A is in the ligand (specified in movableatoms)
882:    IF (STEEREDMINT) THEN831:    IF (STEEREDMINT) THEN
883:       do i=1,nmovableatoms832:       do i=1,nmovableatoms
884:          j=movableatomlist(i)833:          j=movableatomlist(i)
885:          IF (J.EQ.SMINATOMA) GOTO 678834:          IF (J.EQ.SMINATOMA) GOTO 678
1047: CLOSE(20)996: CLOSE(20)
1048: END SUBROUTINE A9DUMPPDB997: END SUBROUTINE A9DUMPPDB
1049: 998: 
1050: SUBROUTINE AMBERFINALIO(nsave,nunit,node,ostring,filth,coords2)999: SUBROUTINE AMBERFINALIO(nsave,nunit,node,ostring,filth,coords2)
1051: ! writes final geometries into AMBER rst and pdb files1000: ! writes final geometries into AMBER rst and pdb files
1052: 1001: 
1053: use modamber91002: use modamber9
1054: 1003: 
1055: implicit none1004: implicit none
1056: 1005: 
1057: integer nsave, i1, nunit,filth,resnumber,currentresidue,node,GETUNIT1006: integer nsave, i1, nunit,filth,resnumber,currentresidue,node
1058: character(len=40) restrt1, restrt2, pdbfile, xyzfile1007: character(len=40) restrt1, restrt2, pdbfile, xyzfile
1059: character(len=10) nsavestring,i1char1008: character(len=10) nsavestring,i1char
1060: character(len=80) ostring,nodestring1009: character(len=80) ostring,nodestring
1061: double precision coords2(3*natom)1010: double precision coords2(3*natom)
1062: 1011: 
1063: ! nunit=25 for GMIN output (lowest.xyz)1012: ! nunit=25 for GMIN output (lowest.xyz)
1064: ! nunit=20 for OPTIM output (points.final.xyz)1013: ! nunit=20 for OPTIM output (points.final.xyz)
1065: 1014: 
1066: ! write(MYUNITNEW,*) 'amberfinalio> ostring = ', ostring1015: ! write(MYUNITNEW,*) 'amberfinalio> ostring = ', ostring
1067: 1016: 
1074: !END IF1023: !END IF
1075: 1024: 
1076: restrt1 = restrt1025: restrt1 = restrt
1077: !write(MYUNITNEW,*) 'NATOM', natom1026: !write(MYUNITNEW,*) 'NATOM', natom
1078: 1027: 
1079:     x(lcrd:lcrd+3*natom-1) = coords2(1:3*natom)1028:     x(lcrd:lcrd+3*natom-1) = coords2(1:3*natom)
1080: IF(nunit==25) THEN      ! dump output files for GMIN1029: IF(nunit==25) THEN      ! dump output files for GMIN
1081:     write(pdbfile,'(A)') 'lowest'//trim(adjustl(nsavestring))//'.pdb'1030:     write(pdbfile,'(A)') 'lowest'//trim(adjustl(nsavestring))//'.pdb'
1082:     write(restrt2,'(A)') 'min'//trim(adjustl(nsavestring))//'.rst'1031:     write(restrt2,'(A)') 'min'//trim(adjustl(nsavestring))//'.rst'
1083: ELSE IF(nunit/=20) THEN ! GMIN with MPI1032: ELSE IF(nunit/=20) THEN ! GMIN with MPI
1084:     write(nodestring,'(I6)') node ! MYNODE in main.F1033:     write(nodestring,'(I6)') node-200 ! ambfinalio_node in main.F (GMIN) is 200+MYNODE
1085:     write(pdbfile,'(A)') 'lowest'//trim(adjustl(nsavestring))//'.'//trim(adjustl(nodestring))//'.pdb'1034:     write(pdbfile,'(A)') 'lowest'//trim(adjustl(nsavestring))//'.'//trim(adjustl(nodestring))//'.pdb'
1086:     write(restrt2,'(A)') 'min'//trim(adjustl(nsavestring))//'.'//trim(adjustl(nodestring))//'.rst'1035:     write(restrt2,'(A)') 'min'//trim(adjustl(nsavestring))//'.'//trim(adjustl(nodestring))//'.rst'
1087: ELSE                    ! nunit=20, dump files for OPTIM1036: ELSE                    ! nunit=20, dump files for OPTIM
1088:  IF(FILTH.EQ.0) THEN1037:  IF(FILTH.EQ.0) THEN
1089:     open(unit=nunit,file='points.final.xyz',status='unknown')1038:     open(unit=nunit,file='points.final.xyz',status='unknown')
1090:     write(nunit,'(I6)') natom1039:     write(nunit,'(I6)') natom
1091:     write(nunit,'(A)') 'points.final'1040:     write(nunit,'(A)') 'points.final'
1092:     write(pdbfile,'(A)') 'points.final.pdb'1041:     write(pdbfile,'(A)') 'points.final.pdb'
1093:     write(restrt2,'(A)') 'points.final.rst'1042:     write(restrt2,'(A)') 'points.final.rst'
1094:  ELSE1043:  ELSE
1095:     open(unit=nunit,file='points.final.xyz.'//trim(adjustl(ostring)),status='unknown')1044:     open(unit=nunit,file='points.final.xyz.'//trim(adjustl(ostring)),status='unknown')
1096:     write(nunit,'(I6)') natom1045:     write(nunit,'(I6)') natom
1097:     write(nunit,'(A)') 'points.final.'//trim(adjustl(ostring))1046:     write(nunit,'(A)') 'points.final.'//trim(adjustl(ostring))
1098:     write(pdbfile,'(A)') 'points.final.pdb.'//trim(adjustl(ostring))1047:     write(pdbfile,'(A)') 'points.final.pdb.'//trim(adjustl(ostring))
1099:     write(restrt2,'(A)') 'points.final.rst.'//trim(adjustl(ostring))1048:     write(restrt2,'(A)') 'points.final.rst.'//trim(adjustl(ostring))
1100:  END IF1049:  END IF
1101: END IF1050: END IF
1102: 1051: 
1103:     restrt = restrt21052:     restrt = restrt2
1104: !WRITE(nunit,*) 'nsave=', nsave, trim(adjustl(pdbfile))1053: 
1105: !WRITE(nunit,*) coords2(1:3) 
1106: !WRITE(nunit,*) x(lcrd:lcrd+2) 
1107: DO i1=1,natom1054: DO i1=1,natom
1108:         WRITE(nunit,'(A4,3F20.10)') ih(m04+i1-1),x(lcrd+3*i1-3),x(lcrd+3*i1-2),x(lcrd+3*i1-1) 1055:         WRITE(nunit,'(A4,3F20.10)') ih(m04+i1-1),x(lcrd+3*i1-3),x(lcrd+3*i1-2),x(lcrd+3*i1-1) 
1109: END DO1056: END DO
1110: ambpdb_unit=GETUNIT()1057: 
1111: open(unit=ambpdb_unit,file=trim(adjustl(pdbfile)),status='unknown')1058: open(unit=ambpdb_unit,file=trim(adjustl(pdbfile)),status='unknown')
1112: 1059: 
1113: ! write out the coordinates to the pdb files1060: ! write out the coordinates to the pdb files
1114: currentresidue=01061: currentresidue=0
1115: DO i1=1,natom1062: DO i1=1,natom
1116: !write(i1char,'(I10)') i11063: !write(i1char,'(I10)') i1
1117: !    resnumber = currentresidue1064: !    resnumber = currentresidue
1118:  IF(i1==1) THEN1065:  IF(i1==1) THEN
1119:     resnumber=11066:     resnumber=1
1120:  END IF1067:  END IF


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0