module MDyn !----------------------------------------------------------------------------- use io_units use names use math use md_calc use geometry_data use io_base use geometry use energy use MD_data use REMD implicit none !----------------------------------------------------------------------------- ! common.MD ! common /mdgrad/ in module.energy ! common /back_constr/ in module.energy ! common /qmeas/ in module.energy ! common /mdpar/ ! common /MDcalc/ ! common /lagrange/ real(kind=8),dimension(:),allocatable :: d_t_work,& d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6) real(kind=8),dimension(:,:),allocatable :: d_t_new,& d_a_old,d_a_short!,d_a !(3,0:MAXRES2) ! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES) ! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,& ! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2) ! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2) ! integer :: dimen,dimen1,dimen3 ! integer :: lang,count_reset_moment,count_reset_vel ! logical :: reset_moment,reset_vel,rattle,RESPA ! common /inertia/ ! common /langevin/ ! real(kind=8) :: rwat,etawat,stdfp,cPoise ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1) ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp) real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES) !----------------------------------------------------------------------------- ! 'sizes.i' ! ! ! ################################################### ! ## COPYRIGHT (C) 1992 by Jay William Ponder ## ! ## All Rights Reserved ## ! ################################################### ! ! ############################################################# ! ## ## ! ## sizes.i -- parameter values to set array dimensions ## ! ## ## ! ############################################################# ! ! ! "sizes.i" sets values for critical array dimensions used ! throughout the software; these parameters will fix the size ! of the largest systems that can be handled; values too large ! for the computer's memory and/or swap space to accomodate ! will result in poor performance or outright failure ! ! parameter: maximum allowed number of: ! ! maxatm atoms in the molecular system ! maxval atoms directly bonded to an atom ! maxgrp ! user-defined groups of atoms ! maxtyp force field atom type definitions ! maxclass force field atom class definitions ! maxkey lines in the keyword file ! maxrot bonds for torsional rotation ! maxvar optimization variables (vector storage) ! maxopt optimization variables (matrix storage) ! maxhess off-diagonal Hessian elements ! maxlight sites for method of lights neighbors ! maxvib vibrational frequencies ! maxgeo distance geometry points ! maxcell unit cells in replicated crystal ! maxring 3-, 4-, or 5-membered rings ! maxfix geometric restraints ! maxbio biopolymer atom definitions ! maxres residues in the macromolecule ! maxamino amino acid residue types ! maxnuc nucleic acid residue types ! maxbnd covalent bonds in molecular system ! maxang bond angles in molecular system ! maxtors torsional angles in molecular system ! maxpi atoms in conjugated pisystem ! maxpib covalent bonds involving pisystem ! maxpit torsional angles involving pisystem ! ! !el integer maxatm,maxval,maxgrp !el integer maxtyp,maxclass,maxkey !el integer maxrot,maxopt !el integer maxhess,maxlight,maxvib !el integer maxgeo,maxcell,maxring !el integer maxfix,maxbio !el integer maxamino,maxnuc,maxbnd !el integer maxang,maxtors,maxpi !el integer maxpib,maxpit ! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres ! integer,parameter :: maxval=8 ! integer,parameter :: maxgrp=1000 ! integer,parameter :: maxtyp=3000 ! integer,parameter :: maxclass=500 ! integer,parameter :: maxkey=10000 ! integer,parameter :: maxrot=1000 ! integer,parameter :: maxopt=1000 ! integer,parameter :: maxhess=1000000 ! integer :: maxlight !=8*maxatm ! integer,parameter :: maxvib=1000 ! integer,parameter :: maxgeo=1000 ! integer,parameter :: maxcell=10000 ! integer,parameter :: maxring=10000 ! integer,parameter :: maxfix=10000 ! integer,parameter :: maxbio=10000 ! integer,parameter :: maxamino=31 ! integer,parameter :: maxnuc=12 ! integer :: maxbnd !=2*maxatm ! integer :: maxang !=3*maxatm ! integer :: maxtors !=4*maxatm ! integer,parameter :: maxpi=100 ! integer,parameter :: maxpib=2*maxpi ! integer,parameter :: maxpit=4*maxpi !----------------------------------------------------------------------------- ! Maximum number of seed ! integer,parameter :: max_seed=1 !----------------------------------------------------------------------------- real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres ! common /stochcalc/ stochforcvec !----------------------------------------------------------------------------- ! common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres !----------------------------------------------------------------------------- ! common /syfek/ subroutines: friction_force,setup_fricmat !el real(kind=8),dimension(:),allocatable :: gamvec !(MAXRES6) or (MAXRES2) !----------------------------------------------------------------------------- ! common /przechowalnia/ subroutines: friction_force,setup_fricmat real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres !----------------------------------------------------------------------------- ! common /przechowalnia/ subroutine: setup_fricmat #ifndef LBFGS real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres) #endif !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains !----------------------------------------------------------------------------- ! brown_step.f !----------------------------------------------------------------------------- subroutine brown_step(itime) !------------------------------------------------ ! Perform a single Euler integration step of Brownian dynamics !------------------------------------------------ ! implicit real*8 (a-h,o-z) use comm_gucio use control, only: tcpu use control_data use energy_data ! use io_conf, only:cartprint ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres integer :: rstcount !ilen, !el external ilen !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres) real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres ! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres) ! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres ! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres !el common /stochcalc/ stochforcvec !el real(kind=8),dimension(3) :: cm !el !el common /gucio/ cm integer :: itime logical :: lprn = .false.,lprn1 = .false. integer :: maxiter = 5 real(kind=8) :: difftol = 1.0d-5 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0 integer :: i,j,nbond,k,ind,ind1,iter integer :: nres2,nres6 logical :: osob nres2=2*nres nres6=6*nres ! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2)) ! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2)) ! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2)) ! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2)) ! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2)) ! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres)) if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres nbond=nct-nnt do i=nnt,nct if (itype(i,1).ne.10) nbond=nbond+1 enddo ! if (lprn1) then write (iout,*) "Generalized inverse of fricmat" call matout(dimen,dimen,nres6,nres6,fricmat) endif do i=1,dimen do j=1,nbond Bmat(i,j)=0.0d0 enddo enddo ind=3 ind1=0 do i=nnt,nct-1 ind1=ind1+1 do j=1,3 Bmat(ind+j,ind1)=dC_norm(j,i) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i,1).ne.10) then ind1=ind1+1 do j=1,3 Bmat(ind+j,ind1)=dC_norm(j,i+nres) enddo ind=ind+3 endif enddo if (lprn1) then write (iout,*) "Matrix Bmat" call MATOUT(nbond,dimen,nres6,nres6,Bmat) endif do i=1,dimen do j=1,nbond GBmat(i,j)=0.0d0 do k=1,dimen GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j) enddo enddo enddo if (lprn1) then write (iout,*) "Matrix GBmat" call MATOUT(nbond,dimen,nres6,nres2,Gbmat) endif do i=1,nbond do j=1,nbond Cmat_(i,j)=0.0d0 do k=1,dimen Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j) enddo enddo enddo if (lprn1) then write (iout,*) "Matrix Cmat" call MATOUT(nbond,nbond,nres2,nres2,Cmat_) endif call matinvert(nbond,nres2,Cmat_,Cinv,osob) if (lprn1) then write (iout,*) "Matrix Cinv" call MATOUT(nbond,nbond,nres2,nres2,Cinv) endif do i=1,dimen do j=1,nbond Tmat(i,j)=0.0d0 do k=1,nbond Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j) enddo enddo enddo if (lprn1) then write (iout,*) "Matrix Tmat" call MATOUT(nbond,dimen,nres6,nres2,Tmat) endif do i=1,dimen do j=1,dimen if (i.eq.j) then Pmat(i,j)=1.0d0 else Pmat(i,j)=0.0d0 endif do k=1,nbond Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k) enddo enddo enddo if (lprn1) then write (iout,*) "Matrix Pmat" call MATOUT(dimen,dimen,nres6,nres6,Pmat) endif do i=1,dimen Td(i)=0.0d0 ind=0 do k=nnt,nct-1 ind=ind+1 Td(i)=Td(i)+vbl*Tmat(i,ind) enddo do k=nnt,nct if (itype(k,1).ne.10) then ind=ind+1 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind) endif enddo enddo if (lprn1) then write (iout,*) "Vector Td" do i=1,dimen write (iout,'(i5,f10.5)') i,Td(i) enddo endif call stochastic_force(stochforcvec) if (lprn) then write (iout,*) "stochforcvec" do i=1,dimen write (iout,*) i,stochforcvec(i) enddo endif do j=1,3 zapas(j)=-gcart(j,0)+stochforcvec(j) d_t_work(j)=d_t(j,0) dC_work(j)=dC_old(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 ind=ind+1 zapas(ind)=-gcart(j,i)+stochforcvec(ind) dC_work(ind)=dC_old(j,i) enddo enddo do i=nnt,nct if (itype(i,1).ne.10) then do j=1,3 ind=ind+1 zapas(ind)=-gxcart(j,i)+stochforcvec(ind) dC_work(ind)=dC_old(j,i+nres) enddo endif enddo if (lprn) then write (iout,*) "Initial d_t_work" do i=1,dimen write (iout,*) i,d_t_work(i) enddo endif do i=1,dimen d_t_work(i)=0.0d0 do j=1,dimen d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j) enddo enddo do i=1,dimen zapas(i)=Td(i) do j=1,dimen zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time) enddo enddo if (lprn1) then write (iout,*) "Final d_t_work and zapas" do i=1,dimen write (iout,*) i,d_t_work(i),zapas(i) enddo endif do j=1,3 d_t(j,0)=d_t_work(j) dc(j,0)=zapas(j) dc_work(j)=dc(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_work(i) dc(j,i)=zapas(ind+j) dc_work(ind+j)=dc(j,i) enddo ind=ind+3 enddo do i=nnt,nct do j=1,3 d_t(j,i+nres)=d_t_work(ind+j) dc(j,i+nres)=zapas(ind+j) dc_work(ind+j)=dc(j,i+nres) enddo ind=ind+3 enddo if (lprn) then call chainbuild_cart write (iout,*) "Before correction for rotational lengthening" write (iout,*) "New coordinates",& " and differences between actual and standard bond lengths" ind=0 do i=nnt,nct-1 ind=ind+1 xx=vbld(i+1)-vbl write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i),j=1,3),xx enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 xx=vbld(i+nres)-vbldsc0(1,itype(i,1)) write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i+nres),j=1,3),xx endif enddo endif ! Second correction (rotational lengthening) ! do iter=1,maxiter diffmax=0.0d0 ind=0 do i=nnt,nct-1 ind=ind+1 blen2 = scalar(dc(1,i),dc(1,i)) ppvec(ind)=2*vbl**2-blen2 diffbond=dabs(vbl-dsqrt(blen2)) if (diffbond.gt.diffmax) diffmax=diffbond if (ppvec(ind).gt.0.0d0) then ppvec(ind)=dsqrt(ppvec(ind)) else ppvec(ind)=0.0d0 endif if (lprn) then write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind) endif enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 blen2 = scalar(dc(1,i+nres),dc(1,i+nres)) ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2 diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2)) if (diffbond.gt.diffmax) diffmax=diffbond if (ppvec(ind).gt.0.0d0) then ppvec(ind)=dsqrt(ppvec(ind)) else ppvec(ind)=0.0d0 endif if (lprn) then write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind) endif endif enddo if (lprn) write (iout,*) "iter",iter," diffmax",diffmax if (diffmax.lt.difftol) goto 10 do i=1,dimen Td(i)=0.0d0 do j=1,nbond Td(i)=Td(i)+ppvec(j)*Tmat(i,j) enddo enddo do i=1,dimen zapas(i)=Td(i) do j=1,dimen zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j) enddo enddo do j=1,3 dc(j,0)=zapas(j) dc_work(j)=zapas(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 dc(j,i)=zapas(ind+j) dc_work(ind+j)=zapas(ind+j) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i,1).ne.10) then do j=1,3 dc(j,i+nres)=zapas(ind+j) dc_work(ind+j)=zapas(ind+j) enddo ind=ind+3 endif enddo ! Building the chain from the newly calculated coordinates call chainbuild_cart if(ntwe.ne.0) then if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "Cartesian and internal coordinates: step 1" call cartprint call intout write (iout,'(a)') "Potential forces" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),& (-gxcart(j,i),j=1,3) enddo write (iout,'(a)') "Stochastic forces" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),& (stochforc(j,i+nres),j=1,3) enddo write (iout,'(a)') "Velocities" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo endif endif if (lprn) then write (iout,*) "After correction for rotational lengthening" write (iout,*) "New coordinates",& " and differences between actual and standard bond lengths" ind=0 do i=nnt,nct-1 ind=ind+1 xx=vbld(i+1)-vbl write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i),j=1,3),xx enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 xx=vbld(i+nres)-vbldsc0(1,itype(i,1)) write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i+nres),j=1,3),xx endif enddo endif ! ENDDO ! write (iout,*) "Too many attempts at correcting the bonds" ! stop 10 continue #ifdef MPI tt0 =MPI_Wtime() #else tt0 = tcpu() #endif ! Calculate energy and forces call zerograd call etotal(potEcomp) potE=potEcomp(0)-potEcomp(51) call cartgrad totT=totT+d_time totTafm=totT ! Calculate the kinetic and total energy and the kinetic temperature call kinetic(EK) #ifdef MPI t_enegrad=t_enegrad+MPI_Wtime()-tt0 #else t_enegrad=t_enegrad+tcpu()-tt0 #endif totE=EK+potE kinetic_T=2.0d0/(dimen*Rb)*EK return end subroutine brown_step !----------------------------------------------------------------------------- ! gauss.f !----------------------------------------------------------------------------- subroutine gauss(RO,AP,MT,M,N,*) ! ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION ! RO IS A SQUARE MATRIX ! THE CALCULATED PRODUCT IS STORED IN AP ! ABNORMAL EXIT IF RO IS SINGULAR ! integer :: MT, M, N, M1,I,J,IM,& I1,MI,MI1 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y integer :: k ! real(kind=8) :: if(M.ne.1)goto 10 X=RO(1,1) if(dabs(X).le.1.0D-13) return 1 X=1.0/X do 16 I=1,N 16 AP(1,I)=AP(1,I)*X return 10 continue M1=M-1 DO 1 I=1,M1 IM=I RM=DABS(RO(I,I)) I1=I+1 do 2 J=I1,M if(DABS(RO(J,I)).LE.RM) goto 2 RM=DABS(RO(J,I)) IM=J 2 continue If(IM.eq.I)goto 17 do 3 J=1,N PR=AP(I,J) AP(I,J)=AP(IM,J) 3 AP(IM,J)=PR do 4 J=I,M PR=RO(I,J) RO(I,J)=RO(IM,J) 4 RO(IM,J)=PR 17 X=RO(I,I) if(dabs(X).le.1.0E-13) return 1 X=1.0/X do 5 J=1,N 5 AP(I,J)=X*AP(I,J) do 6 J=I1,M 6 RO(I,J)=X*RO(I,J) do 7 J=I1,M Y=RO(J,I) do 8 K=1,N 8 AP(J,K)=AP(J,K)-Y*AP(I,K) do 9 K=I1,M 9 RO(J,K)=RO(J,K)-Y*RO(I,K) 7 continue 1 continue X=RO(M,M) if(dabs(X).le.1.0E-13) return 1 X=1.0/X do 11 J=1,N 11 AP(M,J)=X*AP(M,J) do 12 I=1,M1 MI=M-I MI1=MI+1 do 14 J=1,N X=AP(MI,J) do 15 K=MI1,M 15 X=X-AP(K,J)*RO(MI,K) 14 AP(MI,J)=X 12 continue return end subroutine gauss !----------------------------------------------------------------------------- ! kinetic_lesyng.f #ifdef FIVEDIAG subroutine kinetic(KE_total) !c---------------------------------------------------------------- !c This subroutine calculates the total kinetic energy of the chain !c----------------------------------------------------------------- !c 3/5/2020 AL Corrected for multichain systems, no fake peptide groups !c inside, implemented with five-diagonal inertia matrix use energy_data implicit none real(kind=8):: KE_total,KEt_p,KEt_sc,KEr_p,KEr_sc,mag1,mag2 integer i,j,k,iti,mnum real(kind=8),dimension(3) :: incr,v KEt_p=0.0d0 KEt_sc=0.0d0 KEr_p=0.0D0 KEr_sc=0.0D0 !c write (iout,*) "ISC",(isc(itype(i)),i=1,nres) !c The translational part for peptide virtual bonds do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct-1 !czy na pewno nct-1?? mnum=molnum(i) !c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3 !c Skip dummy peptide groups if (itype(i,mnum).ne.ntyp1_molec(mnum)& .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) enddo if (mnum.eq.5) mp(mnum)=0.0d0 ! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) !c write (iout,*) "Kinetic trp:",i,(v(j),j=1,3) vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) endif do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo enddo !c write(iout,*) 'KEt_p', KEt_p !c The translational part for the side chain virtual bond !c Only now we can initialize incr with zeros. It must be equal !c to the velocities of the first Calpha. do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) if (mnum.eq.5) iti=itype(i,mnum) ! if (itype(i,mnum).eq.ntyp1_molec(mnum)) then ! do j=1,3 ! incr(j)=d_t(j,i) ! enddo ! endif if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& .or.mnum.ge.3) then do j=1,3 v(j)=incr(j) enddo else do j=1,3 v(j)=incr(j)+d_t(j,nres+i) enddo endif ! if (mnum.ne.5) then ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3) KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) ! endif do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo enddo ! goto 111 ! write(iout,*) 'KEt_sc', KEt_sc ! The part due to stretching and rotation of the peptide groups do i=nnt,nct-1 mnum=molnum(i) if (itype(i,mnum).ne.ntyp1_molec(mnum)& .and.itype(i+1,mnum).ne.ntyp1_molec(mnum)) then if (mnum.eq.5) Ip(mnum)=0.0 ! write (iout,*) "i",i ! write (iout,*) "i",i," mag1",mag1," mag2",mag2 do j=1,3 incr(j)=d_t(j,i) enddo !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3) KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) & +incr(3)*incr(3)) endif enddo !c goto 111 !c write(iout,*) 'KEr_p', KEr_p !c The rotational part of the side chain virtual bond do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) if (itype(i,1).ne.10.and.itype(i,mnum).ne.ntyp1_molec(mnum)& .and.mnum.lt.3) then do j=1,3 incr(j)=d_t(j,nres+i) enddo ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3) KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+& incr(3)*incr(3)) endif enddo !c The total kinetic energy 111 continue ! write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p, & ! ' KEr_sc', KEr_sc KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc) !c write (iout,*) "KE_total",KE_total return end subroutine kinetic #else !----------------------------------------------------------------------------- subroutine kinetic(KE_total) !---------------------------------------------------------------- ! This subroutine calculates the total kinetic energy of the chain !----------------------------------------------------------------- use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.MD' ! include 'COMMON.IOUNITS' real(kind=8) :: KE_total,mscab integer :: i,j,k,iti,mnum,term real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),& mag1,mag2,v(3) #ifdef DEBUG write (iout,*) "Velocities, kietic" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo #endif KEt_p=0.0d0 KEt_sc=0.0d0 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres) ! The translational part for peptide virtual bonds do j=1,3 incr(j)=d_t(j,0) enddo term=nct-1 ! if (molnum(nct).gt.3) term=nct do i=nnt,term mnum=molnum(i) if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum) if (mnum.gt.4) then do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) enddo else do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) enddo endif vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo enddo ! write(iout,*) 'KEt_p', KEt_p ! The translational part for the side chain virtual bond ! Only now we can initialize incr with zeros. It must be equal ! to the velocities of the first Calpha. do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) ! if (mnum.ge.4) then ! mscab=0.0d0 ! else mscab=msc(iti,mnum) ! endif ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& .or.(mnum.ge.4)) then do j=1,3 v(j)=incr(j) enddo else do j=1,3 v(j)=incr(j)+d_t(j,nres+i) enddo endif ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3) KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo enddo ! goto 111 ! write(iout,*) 'KEt_sc', KEt_sc ! The part due to stretching and rotation of the peptide groups KEr_p=0.0D0 do i=nnt,nct-1 mnum=molnum(i) ! write (iout,*) "i",i ! write (iout,*) "i",i," mag1",mag1," mag2",mag2 do j=1,3 incr(j)=d_t(j,i) enddo ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3) KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) & +incr(3)*incr(3)) enddo ! goto 111 ! write(iout,*) 'KEr_p', KEr_p ! The rotational part of the side chain virtual bond KEr_sc=0.0D0 do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then do j=1,3 incr(j)=d_t(j,nres+i) enddo ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3) KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ & incr(3)*incr(3)) endif enddo ! The total kinetic energy 111 continue ! write(iout,*) 'KEr_sc', KEr_sc KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc) ! write (iout,*) "KE_total",KE_total return end subroutine kinetic !----------------------------------------------------------------------------- #endif subroutine kinetic_CASC(KE_total) !c---------------------------------------------------------------- !c Compute the kinetic energy of the system using the Calpha-SC !c coordinate system !c----------------------------------------------------------------- implicit none double precision KE_total integer i,j,k,iti,ichain,innt,inct,mnum double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),& mag1,mag2,v(3) #ifdef FIVEDIAG KEt_p=0.0d0 KEt_sc=0.0d0 KEr_p=0.0D0 KEr_sc=0.0D0 !c write (iout,*) "ISC",(isc(itype(i)),i=1,nres) !c The translational part for peptide virtual bonds do ichain=1,nchain innt=chain_border(1,ichain) inct=chain_border(2,ichain) !c write (iout,*) "Kinetic_CASC chain",ichain," innt",innt, !c & " inct",inct do i=innt,inct-1 mnum=molnum(i) if (mnum.eq.5) mp(mnum)=0.0d0 ! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) !c write (iout,*) i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3) do j=1,3 v(j)=0.5d0*(d_t(j,i)+d_t(j,i+1)) enddo !c write (iout,*) "Kinetic trp i",i," v",(v(j),j=1,3) KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) enddo !c write(iout,*) 'KEt_p', KEt_p !c The translational part for the side chain virtual bond !c Only now we can initialize incr with zeros. It must be equal !c to the velocities of the first Calpha. do i=innt,inct mnum=molnum(i) if (mnum.eq.5) then iti=itype(i,mnum) else iti=iabs(itype(i,mnum)) endif if (itype(i,1).eq.10.or.mnum.ge.3.or. itype(i,mnum).eq.ntyp1_molec(mnum)) then !c write (iout,*) i,iti,(d_t(j,i),j=1,3) do j=1,3 v(j)=d_t(j,i) enddo else !c write (iout,*) i,iti,(d_t(j,nres+i),j=1,3) do j=1,3 v(j)=d_t(j,nres+i) enddo endif !c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) !c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3) KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) enddo !c goto 111 !c write(iout,*) 'KEt_sc', KEt_sc !c The part due to stretching and rotation of the peptide groups do i=innt,inct-1 mnum=molnum(i) do j=1,3 incr(j)=d_t(j,i+1)-d_t(j,i) enddo if (mnum.eq.5) Ip(mnum)=0.0d0 !c write (iout,*) i,(incr(j),j=1,3) !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3) KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2)& +incr(3)*incr(3)) enddo !c goto 111 !c write(iout,*) 'KEr_p', KEr_p !c The rotational part of the side chain virtual bond do i=innt,inct mnum=molnum(i) ! iti=iabs(itype(i,mnum)) if (mnum.eq.5) then iti=itype(i,mnum) else iti=iabs(itype(i,mnum)) endif ! if (iti.ne.10.and.mnum.lt.3) then if (itype(i,1).ne.10.and.mnum.lt.3.and. itype(i,mnum).ne.ntyp1_molec(mnum)) then do j=1,3 incr(j)=d_t(j,nres+i)-d_t(j,i) enddo !c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3) KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+& incr(3)*incr(3)) endif enddo enddo ! ichain !c The total kinetic energy 111 continue !c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p, !c & ' KEr_sc', KEr_sc KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc) !c write (iout,*) "KE_total",KE_tota #else write (iout,*) "Need to compile with -DFIVEDIAG to use this sub!" stop #endif return end subroutine kinetic_CASC ! MD_A-MTS.F !----------------------------------------------------------------------------- subroutine MD !------------------------------------------------ ! The driver for molecular dynamics subroutines !------------------------------------------------ use comm_gucio ! use MPI use control, only:tcpu,ovrtim ! use io_comm, only:ilen use control_data use compare, only:secondary2,hairpin use io, only:cartout,statout ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI include "mpif.h" integer :: IERROR,ERRCODE #endif ! include 'COMMON.SETUP' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' ! include 'COMMON.HAIRPIN' real(kind=8),dimension(3) :: L,vcm,boxx #ifdef VOUT real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres #endif integer :: rstcount !ilen, !el external ilen character(len=50) :: tytul !el common /gucio/ cm integer :: i,j,nharp integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3) ! logical :: ovrtim real(kind=8) :: tt0,scalfac integer :: nres2,itime nres2=2*nres print *, "ENTER MD" boxx(1)=boxxsize boxx(2)=boxysize boxx(3)=boxzsize ! #ifdef MPI print *,"MY tmpdir",tmpdir,ilen(tmpdir) if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" & //liczba(:ilen(liczba))//'.rst') #else if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst') #endif t_MDsetup=0.0d0 t_langsetup=0.0d0 t_MD=0.0d0 t_enegrad=0.0d0 t_sdsetup=0.0d0 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started" #ifdef MPI tt0=MPI_Wtime() #else tt0 = tcpu() #endif print *,"just befor setup matix",nres ! Determine the inverse of the inertia matrix. call setup_MD_matrices ! Initialize MD print *,"AFTER SETUP MATRICES" call init_MD print *,"AFTER INIT MD" #ifdef MPI t_MDsetup = MPI_Wtime()-tt0 #else t_MDsetup = tcpu()-tt0 #endif rstcount=0 ! Entering the MD loop #ifdef MPI tt0 = MPI_Wtime() #else tt0 = tcpu() #endif if (lang.eq.2 .or. lang.eq.3) then #ifndef LANG0 call setup_fricmat if (lang.eq.2) then call sd_verlet_p_setup else call sd_verlet_ciccotti_setup endif do i=1,dimen do j=1,dimen pfric0_mat(i,j,0)=pfric_mat(i,j) afric0_mat(i,j,0)=afric_mat(i,j) vfric0_mat(i,j,0)=vfric_mat(i,j) prand0_mat(i,j,0)=prand_mat(i,j) vrand0_mat1(i,j,0)=vrand_mat1(i,j) vrand0_mat2(i,j,0)=vrand_mat2(i,j) enddo enddo flag_stoch(0)=.true. do i=1,maxflag_stoch flag_stoch(i)=.false. enddo #else write (iout,*) & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0" #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #endif stop #endif else if (lang.eq.1 .or. lang.eq.4) then print *,"before setup_fricmat" call setup_fricmat print *,"after setup_fricmat" endif #ifdef MPI t_langsetup=MPI_Wtime()-tt0 tt0=MPI_Wtime() #else t_langsetup=tcpu()-tt0 tt0=tcpu() #endif do itime=1,n_timestep if (large) print *,itime,ntwe if (ovrtim()) exit if (large.and. mod(itime,ntwe).eq.0) & write (iout,*) "itime",itime rstcount=rstcount+1 if (lang.gt.0 .and. surfarea .and. & mod(itime,reset_fricmat).eq.0) then if (lang.eq.2 .or. lang.eq.3) then #ifndef LANG0 call setup_fricmat if (lang.eq.2) then call sd_verlet_p_setup else call sd_verlet_ciccotti_setup endif do i=1,dimen do j=1,dimen pfric0_mat(i,j,0)=pfric_mat(i,j) afric0_mat(i,j,0)=afric_mat(i,j) vfric0_mat(i,j,0)=vfric_mat(i,j) prand0_mat(i,j,0)=prand_mat(i,j) vrand0_mat1(i,j,0)=vrand_mat1(i,j) vrand0_mat2(i,j,0)=vrand_mat2(i,j) enddo enddo flag_stoch(0)=.true. do i=1,maxflag_stoch flag_stoch(i)=.false. enddo #endif else if (lang.eq.1 .or. lang.eq.4) then print *,"before setup_fricmat" call setup_fricmat print *,"after setup_fricmat" endif write (iout,'(a,i10)') & "Friction matrix reset based on surface area, itime",itime endif if (reset_vel .and. tbf .and. lang.eq.0 & .and. mod(itime,count_reset_vel).eq.0) then !WARP WATER do i=1,nres if (molnum(i).eq.5) then call to_box(c(1,i),c(2,i),c(3,i)) do j=1,3 if (c(j,i).le.0) c(j,i)=c(j,i)+boxx(j) enddo endif ! write(iout,*) "COORD",c(1,i),c(2,i),c(3,i) enddo call random_vel write(iout,'(a,f20.2)') & "Velocities reset to random values, time",totT do i=0,2*nres do j=1,3 d_t_old(j,i)=d_t(j,i) enddo enddo endif if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then call inertia_tensor call vcm_vel(vcm) do j=1,3 d_t(j,0)=d_t(j,0)-vcm(j) enddo call kinetic(EK) kinetic_T=2.0d0/(dimen3*Rb)*EK scalfac=dsqrt(T_bath/kinetic_T) write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT do i=0,2*nres do j=1,3 d_t_old(j,i)=scalfac*d_t(j,i) enddo enddo endif if (lang.ne.4) then if (RESPA) then ! Time-reversible RESPA algorithm ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992) call RESPA_step(itime) else ! Variable time step algorithm. if (large) print *,"before verlet_step" call velverlet_step(itime) if (large) print *,"after verlet_step" endif else #ifdef BROWN call brown_step(itime) #else print *,"Brown dynamics not here!" #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #endif stop #endif endif itime_mat=itime if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0) then ! call returnbox call statout(itime) ! call returnbox ! call check_ecartint endif #ifdef VOUT do j=1,3 v_work(j)=d_t(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 ind=ind+1 v_work(ind)=d_t(j,i) enddo enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.lt.4) then do j=1,3 ind=ind+1 v_work(ind)=d_t(j,i+nres) enddo endif enddo write (66,'(80f10.5)') & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres) do i=1,ind v_transf(i)=0.0d0 do j=1,ind v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j) enddo v_transf(i)= v_transf(i)*dsqrt(geigen(i)) enddo write (67,'(80f10.5)') (v_transf(i),i=1,ind) #endif endif if (mod(itime,ntwx).eq.0) then call returnbox call enerprint(potEcomp) write (tytul,'("time",f8.2)') totT if(mdpdb) then write(iout,*) "before hairpin" call hairpin(.true.,nharp,iharp) write(iout,*) "before secondary" call secondary2(.true.) write(iout,*) "before pdbout" call pdbout(potE,tytul,ipdb) ! call enerprint(potEcomp) else call cartout(totT) endif if (fodson) then write(iout,*) "starting fodstep" call fodstep(nfodstep) write(iout,*) "after fodstep" call statout(itime) if(mdpdb) then call hairpin(.true.,nharp,iharp) call secondary2(.true.) call pdbout(potE,tytul,ipdb) else call cartout(totT) endif endif endif if (rstcount.eq.1000.or.itime.eq.n_timestep) then open(irest2,file=rest2name,status='unknown') write(irest2,*) totT,EK,potE,totE,t_bath totTafm=totT ! AL 4/17/17: Now writing d_t(0,:) too do i=0,2*nres write (irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo ! AL 4/17/17: Now writing d_c(0,:) too do i=0,2*nres write (irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo close(irest2) rstcount=0 endif enddo #ifdef MPI t_MD=MPI_Wtime()-tt0 #else t_MD=tcpu()-tt0 #endif write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') & ' Timing ',& 'MD calculations setup:',t_MDsetup,& 'Energy & gradient evaluation:',t_enegrad,& 'Stochastic MD setup:',t_langsetup,& 'Stochastic MD step setup:',t_sdsetup,& 'MD steps:',t_MD write (iout,'(/28(1h=),a25,27(1h=))') & ' End of MD calculation ' #ifdef TIMING_ENE write (iout,*) "time for etotal",t_etotal," elong",t_elong,& " eshort",t_eshort write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,& " time_fricmatmult",time_fricmatmult," time_fsample ",& time_fsample #endif return end subroutine MD !----------------------------------------------------------------------------- subroutine velverlet_step(itime) !------------------------------------------------------------------------------- ! Perform a single velocity Verlet step; the time step can be rescaled if ! increments in accelerations exceed the threshold !------------------------------------------------------------------------------- ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' use comm_gucio use control, only:tcpu use control_data use minimm, only:minim_dc #ifdef MPI include 'mpif.h' integer :: ierror,ierrcode real(kind=8) :: errcode #endif ! include 'COMMON.SETUP' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' ! include 'COMMON.MUCA' real(kind=8),dimension(3) :: vcm,incr real(kind=8),dimension(3) :: L integer :: count,rstcount !ilen, !el external ilen character(len=50) :: tytul integer :: maxcount_scale = 30 !el common /gucio/ cm !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime logical :: scalel real(kind=8) :: epdrift,tt0,fac_time ! if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres scalel=.true. icount_scale=0 if (lang.eq.1) then call sddir_precalc if (large) print *,"after sddir_precalc" else if (lang.eq.2 .or. lang.eq.3) then #ifndef LANG0 call stochastic_force(stochforcvec) #else write (iout,*) & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0" #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #endif stop #endif endif itime_scal=0 do while (scalel) icount_scale=icount_scale+1 ! write(iout,*) "icount_scale",icount_scale,scalel if (icount_scale.gt.maxcount_scale) then write (iout,*) & "ERROR: too many attempts at scaling down the time step. ",& "amax=",amax,"epdrift=",epdrift,& "damax=",damax,"edriftmax=",edriftmax,& "d_time=",d_time call flush(iout) #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE) #endif stop endif ! First step of the velocity Verlet algorithm if (lang.eq.2) then #ifndef LANG0 call sd_verlet1 #endif else if (lang.eq.3) then #ifndef LANG0 call sd_verlet1_ciccotti #endif else if (lang.eq.1) then call sddir_verlet1 else call verlet1 endif ! Build the chain from the newly calculated coordinates call chainbuild_cart if (rattle) call rattle1 if (ntwe.ne.0) then if (large) then !.and. mod(itime,ntwe).eq.0) then write (iout,*) "Cartesian and internal coordinates: step 1" call cartprint call intout write (iout,*) "dC" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),& (dc(j,i+nres),j=1,3) enddo write (iout,*) "Accelerations" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),& (d_a(j,i+nres),j=1,3) enddo write (iout,*) "Velocities, step 1" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo endif endif #ifdef MPI tt0 = MPI_Wtime() #else tt0 = tcpu() #endif ! Calculate energy and forces call zerograd call etotal(potEcomp) ! AL 4/17/17: Reduce the steps if NaNs occurred. if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then call enerprint(potEcomp) d_time=d_time/10.0 if (icount_scale.gt.15) then write (iout,*) "Tu jest problem",potEcomp(0),d_time ! call gen_rand_conf(1,*335) ! call minim_dc(potEcomp(0),iretcode,100) ! call zerograd ! call etotal(potEcomp) ! write(iout,*) "needed to repara,",potEcomp endif cycle ! 335 write(iout,*) "Failed genrand" ! cycle endif ! end change if (large.and. mod(itime,ntwe).eq.0) & call enerprint(potEcomp) #ifdef TIMING_ENE #ifdef MPI t_etotal=t_etotal+MPI_Wtime()-tt0 #else t_etotal=t_etotal+tcpu()-tt0 #endif #endif potE=potEcomp(0)-potEcomp(51) call cartgrad ! Get the new accelerations call lagrangian #ifdef MPI t_enegrad=t_enegrad+MPI_Wtime()-tt0 #else t_enegrad=t_enegrad+tcpu()-tt0 #endif ! Determine maximum acceleration and scale down the timestep if needed call max_accel amax=amax/(itime_scal+1)**2 call predict_edrift(epdrift) ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1) scalel=.false. ! write (iout,*) "before enter if",scalel,icount_scale if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then ! write(iout,*) "I enter if" ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step scalel=.true. ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) & /dlog(2.0d0)+1 itime_scal=itime_scal+ifac_time ! fac_time=dmin1(damax/amax,0.5d0) fac_time=0.5d0**ifac_time d_time=d_time*fac_time if (lang.eq.2 .or. lang.eq.3) then #ifndef LANG0 ! write (iout,*) "Calling sd_verlet_setup: 1" ! Rescale the stochastic forces and recalculate or restore ! the matrices of tinker integrator if (itime_scal.gt.maxflag_stoch) then if (large) write (iout,'(a,i5,a)') & "Calculate matrices for stochastic step;",& " itime_scal ",itime_scal if (lang.eq.2) then call sd_verlet_p_setup else call sd_verlet_ciccotti_setup endif write (iout,'(2a,i3,a,i3,1h.)') & "Warning: cannot store matrices for stochastic",& " integration because the index",itime_scal,& " is greater than",maxflag_stoch write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",& " integration Langevin algorithm for better efficiency." else if (flag_stoch(itime_scal)) then if (large) write (iout,'(a,i5,a,l1)') & "Restore matrices for stochastic step; itime_scal ",& itime_scal," flag ",flag_stoch(itime_scal) do i=1,dimen do j=1,dimen pfric_mat(i,j)=pfric0_mat(i,j,itime_scal) afric_mat(i,j)=afric0_mat(i,j,itime_scal) vfric_mat(i,j)=vfric0_mat(i,j,itime_scal) prand_mat(i,j)=prand0_mat(i,j,itime_scal) vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal) vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal) enddo enddo else if (large) write (iout,'(2a,i5,a,l1)') & "Calculate & store matrices for stochastic step;",& " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal) if (lang.eq.2) then call sd_verlet_p_setup else call sd_verlet_ciccotti_setup endif flag_stoch(ifac_time)=.true. do i=1,dimen do j=1,dimen pfric0_mat(i,j,itime_scal)=pfric_mat(i,j) afric0_mat(i,j,itime_scal)=afric_mat(i,j) vfric0_mat(i,j,itime_scal)=vfric_mat(i,j) prand0_mat(i,j,itime_scal)=prand_mat(i,j) vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j) vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j) enddo enddo endif fac_time=1.0d0/dsqrt(fac_time) do i=1,dimen stochforcvec(i)=fac_time*stochforcvec(i) enddo #endif else if (lang.eq.1) then ! Rescale the accelerations due to stochastic forces fac_time=1.0d0/dsqrt(fac_time) do i=1,dimen d_as_work(i)=d_as_work(i)*fac_time enddo endif if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') & "itime",itime," Timestep scaled down to ",& d_time," ifac_time",ifac_time," itime_scal",itime_scal else ! Second step of the velocity Verlet algorithm if (lang.eq.2) then #ifndef LANG0 call sd_verlet2 #endif else if (lang.eq.3) then #ifndef LANG0 call sd_verlet2_ciccotti #endif else if (lang.eq.1) then call sddir_verlet2 else call verlet2 endif if (rattle) call rattle2 totT=totT+d_time totTafm=totT if (d_time.ne.d_time0) then d_time=d_time0 #ifndef LANG0 if (lang.eq.2 .or. lang.eq.3) then if (large) write (iout,'(a)') & "Restore original matrices for stochastic step" ! write (iout,*) "Calling sd_verlet_setup: 2" ! Restore the matrices of tinker integrator if the time step has been restored do i=1,dimen do j=1,dimen pfric_mat(i,j)=pfric0_mat(i,j,0) afric_mat(i,j)=afric0_mat(i,j,0) vfric_mat(i,j)=vfric0_mat(i,j,0) prand_mat(i,j)=prand0_mat(i,j,0) vrand_mat1(i,j)=vrand0_mat1(i,j,0) vrand_mat2(i,j)=vrand0_mat2(i,j,0) enddo enddo endif #endif endif endif enddo ! Calculate the kinetic and the total energy and the kinetic temperature call kinetic(EK) totE=EK+potE ! diagnostics ! call kinetic1(EK1) ! write (iout,*) "step",itime," EK",EK," EK1",EK1 ! end diagnostics ! Couple the system to Berendsen bath if needed if (tbf .and. lang.eq.0) then call verlet_bath endif kinetic_T=2.0d0/(dimen3*Rb)*EK ! Backup the coordinates, velocities, and accelerations do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0 .and. large) then write (iout,*) "Velocities, step 2" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo endif endif return end subroutine velverlet_step !----------------------------------------------------------------------------- subroutine RESPA_step(itime) !------------------------------------------------------------------------------- ! Perform a single RESPA step. !------------------------------------------------------------------------------- ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' use comm_gucio use comm_cipiszcze ! use MPI use control, only:tcpu use control_data ! use io_conf, only:cartprint #ifdef MPI include 'mpif.h' integer :: IERROR,ERRCODE #endif ! include 'COMMON.SETUP' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' real(kind=8),dimension(0:n_ene) :: energia_short,energia_long real(kind=8),dimension(3) :: L,vcm,incr real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres logical :: PRINT_AMTS_MSG = .false. integer :: count,rstcount !ilen, !el external ilen character(len=50) :: tytul integer :: maxcount_scale = 10 !el common /gucio/ cm !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec integer :: itt,i,j,itsplit,itime logical :: scale !el common /cipiszcze/ itt real(kind=8) :: epdrift,tt0,epdriftmax itt = itt_comm if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres itt=itime if (ntwe.ne.0) then if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "***************** RESPA itime",itime write (iout,*) "Cartesian and internal coordinates: step 0" ! call cartprint call pdbout(0.0d0,"cipiszcze",iout) call intout write (iout,*) "Accelerations from long-range forces" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),& (d_a(j,i+nres),j=1,3) enddo write (iout,*) "Velocities, step 0" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo endif endif ! ! Perform the initial RESPA step (increment velocities) ! write (iout,*) "*********************** RESPA ini" call RESPA_vel if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0 .and. large) then write (iout,*) "Velocities, end" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo endif endif ! Compute the short-range forces #ifdef MPI tt0 =MPI_Wtime() #else tt0 = tcpu() #endif ! 7/2/2009 commented out ! call zerograd ! call etotal_short(energia_short) ! call cartgrad ! call lagrangian ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step do i=0,2*nres do j=1,3 d_a(j,i)=d_a_short(j,i) enddo enddo if (ntwe.ne.0) then if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "energia_short",energia_short(0) write (iout,*) "Accelerations from short-range forces" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),& (d_a(j,i+nres),j=1,3) enddo endif endif #ifdef MPI t_enegrad=t_enegrad+MPI_Wtime()-tt0 #else t_enegrad=t_enegrad+tcpu()-tt0 #endif do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo ! 6/30/08 A-MTS: attempt at increasing the split number do i=0,2*nres do j=1,3 dc_old0(j,i)=dc_old(j,i) d_t_old0(j,i)=d_t_old(j,i) d_a_old0(j,i)=d_a_old(j,i) enddo enddo if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0 ! scale=.true. d_time0=d_time do while (scale) scale=.false. ! write (iout,*) "itime",itime," ntime_split",ntime_split ! Split the time step d_time=d_time0/ntime_split ! Perform the short-range RESPA steps (velocity Verlet increments of ! positions and velocities using short-range forces) ! write (iout,*) "*********************** RESPA split" do itsplit=1,ntime_split if (lang.eq.1) then call sddir_precalc else if (lang.eq.2 .or. lang.eq.3) then #ifndef LANG0 call stochastic_force(stochforcvec) #else write (iout,*) & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0" #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #endif stop #endif endif ! First step of the velocity Verlet algorithm if (lang.eq.2) then #ifndef LANG0 call sd_verlet1 #endif else if (lang.eq.3) then #ifndef LANG0 call sd_verlet1_ciccotti #endif else if (lang.eq.1) then call sddir_verlet1 else call verlet1 endif ! Build the chain from the newly calculated coordinates call chainbuild_cart if (rattle) call rattle1 if (ntwe.ne.0) then if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "***** ITSPLIT",itsplit write (iout,*) "Cartesian and internal coordinates: step 1" call pdbout(0.0d0,"cipiszcze",iout) ! call cartprint call intout write (iout,*) "Velocities, step 1" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo endif endif #ifdef MPI tt0 = MPI_Wtime() #else tt0 = tcpu() #endif ! Calculate energy and forces call zerograd call etotal_short(energia_short) ! AL 4/17/17: Exit itime_split loop when energy goes infinite if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then if (PRINT_AMTS_MSG) & write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split ntime_split=ntime_split*2 if (ntime_split.gt.maxtime_split) then #ifdef MPI write (iout,*) & "Cannot rescue the run; aborting job. Retry with a smaller time step" call flush(iout) call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #else write (iout,*) & "Cannot rescue the run; terminating. Retry with a smaller time step" #endif endif exit endif ! End change if (large.and. mod(itime,ntwe).eq.0) & call enerprint(energia_short) #ifdef TIMING_ENE #ifdef MPI t_eshort=t_eshort+MPI_Wtime()-tt0 #else t_eshort=t_eshort+tcpu()-tt0 #endif #endif call cartgrad ! Get the new accelerations call lagrangian ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array do i=0,2*nres do j=1,3 d_a_short(j,i)=d_a(j,i) enddo enddo if (ntwe.ne.0) then if (large.and. mod(itime,ntwe).eq.0) then write (iout,*)"energia_short",energia_short(0) write (iout,*) "Accelerations from short-range forces" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),& (d_a(j,i+nres),j=1,3) enddo endif endif ! 6/30/08 A-MTS ! Determine maximum acceleration and scale down the timestep if needed call max_accel amax=amax/ntime_split**2 call predict_edrift(epdrift) if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) & write (iout,*) "amax",amax," damax",damax,& " epdrift",epdrift," epdriftmax",epdriftmax ! Exit loop and try with increased split number if the change of ! acceleration is too big if (amax.gt.damax .or. epdrift.gt.edriftmax) then if (ntime_split.lt.maxtime_split) then scale=.true. ntime_split=ntime_split*2 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big exit do i=0,2*nres do j=1,3 dc_old(j,i)=dc_old0(j,i) d_t_old(j,i)=d_t_old0(j,i) d_a_old(j,i)=d_a_old0(j,i) enddo enddo if (PRINT_AMTS_MSG) then write (iout,*) "acceleration/energy drift too large",amax,& epdrift," split increased to ",ntime_split," itime",itime,& " itsplit",itsplit endif exit else write (iout,*) & "Uh-hu. Bumpy landscape. Maximum splitting number",& maxtime_split,& " already reached!!! Trying to carry on!" endif endif #ifdef MPI t_enegrad=t_enegrad+MPI_Wtime()-tt0 #else t_enegrad=t_enegrad+tcpu()-tt0 #endif ! Second step of the velocity Verlet algorithm if (lang.eq.2) then #ifndef LANG0 call sd_verlet2 #endif else if (lang.eq.3) then #ifndef LANG0 call sd_verlet2_ciccotti #endif else if (lang.eq.1) then call sddir_verlet2 else call verlet2 endif if (rattle) call rattle2 ! Backup the coordinates, velocities, and accelerations do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo enddo enddo ! while scale ! Restore the time step d_time=d_time0 ! Compute long-range forces #ifdef MPI tt0 =MPI_Wtime() #else tt0 = tcpu() #endif call zerograd call etotal_long(energia_long) if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then #ifdef MPI write (iout,*) & "Infinitied/NaNs in energia_long, Aborting MPI job." call flush(iout) call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #else write (iout,*) "Infinitied/NaNs in energia_long, terminating." stop #endif endif if (large.and. mod(itime,ntwe).eq.0) & call enerprint(energia_long) #ifdef TIMING_ENE #ifdef MPI t_elong=t_elong+MPI_Wtime()-tt0 #else t_elong=t_elong+tcpu()-tt0 #endif #endif potE=potEcomp(0)-potEcomp(51) call cartgrad call lagrangian #ifdef MPI t_enegrad=t_enegrad+MPI_Wtime()-tt0 #else t_enegrad=t_enegrad+tcpu()-tt0 #endif ! Compute accelerations from long-range forces if (ntwe.ne.0) then if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "energia_long",energia_long(0) write (iout,*) "Cartesian and internal coordinates: step 2" ! call cartprint call pdbout(0.0d0,"cipiszcze",iout) call intout write (iout,*) "Accelerations from long-range forces" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),& (d_a(j,i+nres),j=1,3) enddo write (iout,*) "Velocities, step 2" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo endif endif ! Compute the final RESPA step (increment velocities) ! write (iout,*) "*********************** RESPA fin" call RESPA_vel ! Compute the complete potential energy do i=0,n_ene potEcomp(i)=energia_short(i)+energia_long(i) enddo potE=potEcomp(0)-potEcomp(51) ! potE=energia_short(0)+energia_long(0) totT=totT+d_time totTafm=totT ! Calculate the kinetic and the total energy and the kinetic temperature call kinetic(EK) totE=EK+potE ! Couple the system to Berendsen bath if needed if (tbf .and. lang.eq.0) then call verlet_bath endif kinetic_T=2.0d0/(dimen3*Rb)*EK ! Backup the coordinates, velocities, and accelerations if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0 .and. large) then write (iout,*) "Velocities, end" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo endif endif return end subroutine RESPA_step !----------------------------------------------------------------------------- subroutine RESPA_vel ! First and last RESPA step (incrementing velocities using long-range ! forces). use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' integer :: i,j,inres,mnum do j=1,3 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time enddo do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time enddo enddo do i=nnt,nct mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time enddo endif enddo return end subroutine RESPA_vel !----------------------------------------------------------------------------- subroutine verlet1 ! Applying velocity Verlet algorithm - step 1 to coordinates use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' real(kind=8) :: adt,adt2 integer :: i,j,inres,mnum #ifdef DEBUG write (iout,*) "VELVERLET1 START: DC" do i=0,nres write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),& (dc(j,i+nres),j=1,3) enddo #endif do j=1,3 adt=d_a_old(j,0)*d_time adt2=0.5d0*adt dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time d_t_new(j,0)=d_t_old(j,0)+adt2 d_t(j,0)=d_t_old(j,0)+adt enddo do i=nnt,nct-1 do j=1,3 adt=d_a_old(j,i)*d_time adt2=0.5d0*adt dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time d_t_new(j,i)=d_t_old(j,i)+adt2 d_t(j,i)=d_t_old(j,i)+adt enddo enddo do i=nnt,nct mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then inres=i+nres do j=1,3 adt=d_a_old(j,inres)*d_time adt2=0.5d0*adt dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time d_t_new(j,inres)=d_t_old(j,inres)+adt2 d_t(j,inres)=d_t_old(j,inres)+adt enddo endif enddo #ifdef DEBUG write (iout,*) "VELVERLET1 END: DC" do i=0,nres write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),& (dc(j,i+nres),j=1,3) enddo #endif return end subroutine verlet1 !----------------------------------------------------------------------------- subroutine verlet2 ! Step 2 of the velocity Verlet algorithm: update velocities use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' integer :: i,j,inres,mnum do j=1,3 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time enddo do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time enddo enddo do i=nnt,nct mnum=molnum(i) ! iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time enddo endif enddo return end subroutine verlet2 !----------------------------------------------------------------------------- subroutine sddir_precalc ! Applying velocity Verlet algorithm - step 1 to coordinates ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' use MPI_data use control_data #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec real(kind=8) :: time00 integer :: i ! ! Compute friction and stochastic forces ! #ifdef MPI time00=MPI_Wtime() if (large) print *,"before friction_force" call friction_force if (large) print *,"after friction_force" time_fric=time_fric+MPI_Wtime()-time00 time00=MPI_Wtime() call stochastic_force(stochforcvec) time_stoch=time_stoch+MPI_Wtime()-time00 #endif ! ! Compute the acceleration due to friction forces (d_af_work) and stochastic ! forces (d_as_work) ! ! call ginv_mult(fric_work, d_af_work) ! call ginv_mult(stochforcvec, d_as_work) #ifdef FIVEDIAG write(iout,*) "forces before fivediaginv" do i=1,dimen*3 write(iout,*) "fricwork",i,fric_work(i) enddo call fivediaginv_mult(dimen,fric_work, d_af_work) call fivediaginv_mult(dimen,stochforcvec, d_as_work) if (large) then write(iout,*),"dimen",dimen do i=1,dimen write(iout,*) "fricwork",fric_work(i), d_af_work(i) write(iout,*) "stochforcevec", stochforcvec(i), d_as_work(i) enddo endif #else call ginv_mult(fric_work, d_af_work) call ginv_mult(stochforcvec, d_as_work) #endif return end subroutine sddir_precalc !----------------------------------------------------------------------------- subroutine sddir_verlet1 ! Applying velocity Verlet algorithm - step 1 to velocities ! use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! Revised 3/31/05 AL: correlation between random contributions to ! position and velocity increments included. real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3) real(kind=8) :: adt,adt2 integer :: i,j,ind,inres,mnum ! ! Add the contribution from BOTH friction and stochastic force to the ! coordinates, but ONLY the contribution from the friction forces to velocities ! do j=1,3 adt=(d_a_old(j,0)+d_af_work(j))*d_time adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,0),d_af_work(j),d_time dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt d_t(j,0)=d_t_old(j,0)+adt enddo ind=3 do i=nnt,nct-1 do j=1,3 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,i),d_af_work(ind+j) dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt d_t(j,i)=d_t_old(j,i)+adt enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) ! iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then inres=i+nres do j=1,3 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,inres),d_af_work(ind+j) dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt d_t(j,inres)=d_t_old(j,inres)+adt enddo ind=ind+3 endif enddo return end subroutine sddir_verlet1 !----------------------------------------------------------------------------- subroutine sddir_verlet2 ! Calculating the adjusted velocities for accelerations ! use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' real(kind=8),dimension(:),allocatable :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0 integer :: i,j,ind,inres,mnum if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) if (.not.allocated(d_as_work1)) allocate(d_as_work1(6*nres)) ! Revised 3/31/05 AL: correlation between random contributions to ! position and velocity increments included. ! The correlation coefficients are calculated at low-friction limit. ! Also, friction forces are now not calculated with new velocities. ! call friction_force call stochastic_force(stochforcvec) ! ! Compute the acceleration due to friction forces (d_af_work) and stochastic ! forces (d_as_work) ! #ifdef FIVEDIAG call fivediaginv_mult(6*nres,stochforcvec, d_as_work1) #else call ginv_mult(stochforcvec, d_as_work1) #endif ! ! Update velocities ! do j=1,3 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time enddo ind=3 do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) ! iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) & +d_af_work(ind+j))+sin60*d_as_work(ind+j) & +cos60*d_as_work1(ind+j))*d_time enddo ind=ind+3 endif enddo return end subroutine sddir_verlet2 !----------------------------------------------------------------------------- subroutine max_accel ! ! Find the maximum difference in the accelerations of the the sites ! at the beginning and the end of the time step. ! use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' real(kind=8),dimension(3) :: aux,accel,accel_old real(kind=8) :: dacc integer :: i,j,mnum do j=1,3 ! aux(j)=d_a(j,0)-d_a_old(j,0) accel_old(j)=d_a_old(j,0) accel(j)=d_a(j,0) enddo amax=0.0d0 do i=nnt,nct ! Backbone if (i.lt.nct) then ! 7/3/08 changed to asymmetric difference do j=1,3 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i)) accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i) accel(j)=accel(j)+0.5d0*d_a(j,i) ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) if (dabs(accel(j)).gt.dabs(accel_old(j))) then dacc=dabs(accel(j)-accel_old(j)) ! write (iout,*) i,dacc if (dacc.gt.amax) amax=dacc endif enddo endif enddo ! Side chains do j=1,3 ! accel(j)=aux(j) accel_old(j)=d_a_old(j,0) accel(j)=d_a(j,0) enddo if (nnt.eq.2) then do j=1,3 accel_old(j)=accel_old(j)+d_a_old(j,1) accel(j)=accel(j)+d_a(j,1) enddo endif do i=nnt,nct mnum=molnum(i) ! iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then do j=1,3 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres) accel_old(j)=accel_old(j)+d_a_old(j,i+nres) accel(j)=accel(j)+d_a(j,i+nres) enddo endif do j=1,3 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) if (dabs(accel(j)).gt.dabs(accel_old(j))) then dacc=dabs(accel(j)-accel_old(j)) ! write (iout,*) "side-chain",i,dacc if (dacc.gt.amax) amax=dacc endif enddo do j=1,3 accel_old(j)=accel_old(j)+d_a_old(j,i) accel(j)=accel(j)+d_a(j,i) ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i) enddo enddo return end subroutine max_accel !----------------------------------------------------------------------------- subroutine predict_edrift(epdrift) ! ! Predict the drift of the potential energy ! use energy_data use control_data, only: lmuca ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.MUCA' real(kind=8) :: epdrift,epdriftij integer :: i,j ! Drift of the potential energy epdrift=0.0d0 do i=nnt,nct ! Backbone if (i.lt.nct) then do j=1,3 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i)) if (lmuca) epdriftij=epdriftij*factor ! write (iout,*) "back",i,j,epdriftij if (epdriftij.gt.epdrift) epdrift=epdriftij enddo endif ! Side chains if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.& molnum(i).lt.4) then do j=1,3 epdriftij= & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i)) if (lmuca) epdriftij=epdriftij*factor ! write (iout,*) "side",i,j,epdriftij if (epdriftij.gt.epdrift) epdrift=epdriftij enddo endif enddo epdrift=0.5d0*epdrift*d_time*d_time ! write (iout,*) "epdrift",epdrift return end subroutine predict_edrift !----------------------------------------------------------------------------- subroutine verlet_bath ! ! Coupling to the thermostat by using the Berendsen algorithm ! use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' real(kind=8) :: T_half,fact integer :: i,j,inres,mnum ! T_half=2.0d0/(dimen3*Rb)*EK fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0)) ! write(iout,*) "T_half", T_half ! write(iout,*) "EK", EK ! write(iout,*) "fact", fact do j=1,3 d_t(j,0)=fact*d_t(j,0) enddo do i=nnt,nct-1 do j=1,3 d_t(j,i)=fact*d_t(j,i) enddo enddo do i=nnt,nct mnum=molnum(i) ! iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then inres=i+nres do j=1,3 d_t(j,inres)=fact*d_t(j,inres) enddo endif enddo return end subroutine verlet_bath !----------------------------------------------------------------------------- subroutine init_MD ! Set up the initial conditions of a MD simulation use comm_gucio use energy_data use control, only:tcpu !el use io_basic, only:ilen use control_data use MPI_data use minimm, only:minim_dc,minimize,sc_move use io_config, only:readrst use io, only:statout use random, only: iran_num ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MP include 'mpif.h' character(len=16) :: form integer :: IERROR,ERRCODE #endif integer :: iranmin,itrial,itmp,n_model_try,k, & i_model integer, dimension(:),allocatable :: list_model_try integer, dimension(0:nodes-1) :: i_start_models ! include 'COMMON.SETUP' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.REMD' real(kind=8),dimension(0:n_ene) :: energia_long,energia_short,energia real(kind=8),dimension(3) :: vcm,incr,L real(kind=8) :: xv,sigv,lowb,highb real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) character(len=256) :: qstr !el integer ilen !el external ilen character(len=50) :: tytul logical :: file_exist !el common /gucio/ cm integer :: i,j,ipos,iq,iw,nft_sc,iretcode,ierr,mnum,itime #ifndef LBFGS integer :: nfun #endif real(kind=8) :: etot,tt0 logical :: fail d_time0=d_time ! write(iout,*) "d_time", d_time ! Compute the standard deviations of stochastic forces for Langevin dynamics ! if the friction coefficients do not depend on surface area if (lang.gt.0 .and. .not.surfarea) then do i=nnt,nct-1 mnum=molnum(i) stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum)) enddo do i=nnt,nct mnum=molnum(i) stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) & *dsqrt(gamsc(iabs(itype(i,mnum)),mnum)) enddo endif ! Open the pdb file for snapshotshots #ifdef MPI if(mdpdb) then if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// & liczba(:ilen(liczba))//".pdb") open(ipdb,& file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) & //".pdb") else #ifdef NOXDR if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// & liczba(:ilen(liczba))//".x") cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) & //".x" #else if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// & liczba(:ilen(liczba))//".cx") cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) & //".cx" #endif endif #else if(mdpdb) then if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb") open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb") else if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx") cartname=prefix(:ilen(prefix))//"_MD.cx" endif #endif if (usampl) then write (qstr,'(256(1h ))') ipos=1 do i=1,nfrag iq = qinfrag(i,iset)*10 iw = wfrag(i,iset)/100 if (iw.gt.0) then if(me.eq.king.or..not.out1file) & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw ipos=ipos+7 endif enddo do i=1,npair iq = qinpair(i,iset)*10 iw = wpair(i,iset)/100 if (iw.gt.0) then if(me.eq.king.or..not.out1file) & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw ipos=ipos+7 endif enddo ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb' #ifdef NOXDR ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x' #else ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx' #endif ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat' endif icg=1 if (rest) then if (restart1file) then if (me.eq.king) & inquire(file=mremd_rst_name,exist=file_exist) write (*,*) me," Before broadcast: file_exist",file_exist #ifdef MPI !el call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,& IERR) #endif !el write (*,*) me," After broadcast: file_exist",file_exist ! inquire(file=mremd_rst_name,exist=file_exist) if(me.eq.king.or..not.out1file) & write(iout,*) "Initial state read by master and distributed" else if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' & //liczba(:ilen(liczba))//'.rst') inquire(file=rest2name,exist=file_exist) endif if(file_exist) then if(.not.restart1file) then if(me.eq.king.or..not.out1file) & write(iout,*) "Initial state will be read from file ",& rest2name(:ilen(rest2name)) call readrst endif call rescale_weights(t_bath) else if(me.eq.king.or..not.out1file)then if (restart1file) then write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),& " does not exist" else write(iout,*) "File ",rest2name(:ilen(rest2name)),& " does not exist" endif write(iout,*) "Initial velocities randomly generated" endif call random_vel totT=0.0d0 totTafm=totT endif else ! Generate initial velocities if(me.eq.king.or..not.out1file) & write(iout,*) "Initial velocities randomly generated" call random_vel totT=0.0d0 totTafm=totT endif ! rest2name = prefix(:ilen(prefix))//'.rst' if(me.eq.king.or..not.out1file)then write (iout,*) "Initial velocities" do i=0,nres write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo ! Zeroing the total angular momentum of the system write(iout,*) "Calling the zero-angular momentum subroutine" endif call inertia_tensor ! Getting the potential energy and forces and velocities and accelerations call vcm_vel(vcm) write (iout,*) "velocity of the center of the mass:" write (iout,*) (vcm(j),j=1,3) do j=1,3 d_t(j,0)=d_t(j,0)-vcm(j) enddo ! Removing the velocity of the center of mass call vcm_vel(vcm) if(me.eq.king.or..not.out1file)then write (iout,*) "vcm right after adjustment:" write (iout,*) (vcm(j),j=1,3) endif ! call chainbuild if ((.not.rest).or.(forceminim)) then if (forceminim) call chainbuild_cart 122 continue if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then if (overlapsc) then print *, 'Calling OVERLAP_SC' call overlap_sc(fail) print *,'after OVERLAP' endif if (searchsc) then print *,'call SC_MOVE' call sc_move(2,nres-1,10,1d10,nft_sc,etot) print *,'SC_move',nft_sc,etot if(me.eq.king.or..not.out1file) & write(iout,*) 'SC_move',nft_sc,etot endif if(dccart)then print *, 'Calling MINIM_DC' call minim_dc(etot,iretcode,nfun) else call geom_to_var(nvar,varia) print *,'Calling MINIMIZE.' call minimize(etot,varia,iretcode,nfun) call var_to_geom(nvar,varia) endif write(iout,*) "just before minimin" call cartprint if(me.eq.king.or..not.out1file) & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun endif write(iout,*) "just after minimin" call cartprint if(iranconf.ne.0) then !c 8/22/17 AL Loop to produce a low-energy random conformation DO iranmin=1,40 if (overlapsc) then if(me.eq.king.or..not.out1file) & write (iout,*) 'Calling OVERLAP_SC' call overlap_sc(fail) endif !endif overlap if (searchsc) then call sc_move(2,nres-1,10,1d10,nft_sc,etot) print *,'SC_move',nft_sc,etot if(me.eq.king.or..not.out1file) & write(iout,*) 'SC_move',nft_sc,etot endif if(dccart)then print *, 'Calling MINIM_DC' call minim_dc(etot,iretcode,nfun) call int_from_cart1(.false.) else call geom_to_var(nvar,varia) print *,'Calling MINIMIZE.' call minimize(etot,varia,iretcode,nfun) call var_to_geom(nvar,varia) endif if(me.eq.king.or..not.out1file) & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun write(iout,*) "just after minimin" call cartprint if (isnan(etot) .or. etot.gt.4.0d6) then write (iout,*) "Energy too large",etot, & " trying another random conformation" do itrial=1,100 itmp=1 call gen_rand_conf(itmp,*30) goto 40 30 write (iout,*) 'Failed to generate random conformation', & ', itrial=',itrial write (*,*) 'Processor:',me, & ' Failed to generate random conformation',& ' itrial=',itrial call intout #ifdef AIX call flush_(iout) #else call flush(iout) #endif enddo write (iout,'(a,i3,a)') 'Processor:',me, & ' error in generating random conformation.' write (*,'(a,i3,a)') 'Processor:',me, & ' error in generating random conformation.' call flush(iout) #ifdef MPI ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode) call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #else stop #endif 40 continue else goto 44 endif ENDDO write (iout,'(a,i3,a)') 'Processor:',me, & ' failed to generate a low-energy random conformation.' write (*,'(a,i3,a,f10.3)') 'Processor:',me, & ' failed to generate a low-energy random conformation.',etot call flush(iout) call intout #ifdef MPI ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode) call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #else stop #endif 44 continue else if (preminim) then if (start_from_model) then n_model_try=0 fail=.true. list_model_try=0 do while (fail .and. n_model_try.lt.nmodel_start) write (iout,*) "n_model_try",n_model_try do i_model=iran_num(1,nmodel_start) do k=1,n_model_try if (i_model.eq.list_model_try(k)) exit enddo if (k.gt.n_model_try) exit enddo n_model_try=n_model_try+1 list_model_try(n_model_try)=i_model if (me.eq.king .or. .not. out1file) & write (iout,*) 'Trying to start from model ',& pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model))) do i=1,2*nres do j=1,3 c(j,i)=chomo(j,i,i_model) enddo enddo call int_from_cart(.true.,.false.) call sc_loc_geom(.false.) dc(:,0)=c(:,1) do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) enddo enddo do i=2,nres-1 do j=1,3 dc(j,i+nres)=c(j,i+nres)-c(j,i) dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) enddo enddo if (me.eq.king.or..not.out1file) then write (iout,*) "Energies before removing overlaps" call etotal(energia(0)) call enerprint(energia(0)) endif ! Remove SC overlaps if requested if (overlapsc) then write (iout,*) 'Calling OVERLAP_SC' call overlap_sc(fail) if (fail) then write (iout,*)& "Failed to remove overlap from model",i_model cycle endif endif if (me.eq.king.or..not.out1file) then write (iout,*) "Energies after removing overlaps" call etotal(energia(0)) call enerprint(energia(0)) endif #ifdef SEARCHSC ! Search for better SC rotamers if requested if (searchsc) then call sc_move(2,nres-1,10,1d10,nft_sc,etot) print *,'SC_move',nft_sc,etot if (me.eq.king.or..not.out1file)& write(iout,*) 'SC_move',nft_sc,etot endif call etotal(energia(0)) #endif enddo call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),& 1,MPI_INTEGER,king,CG_COMM,IERROR) if (n_model_try.gt.nmodel_start .and.& (me.eq.king .or. out1file)) then write (iout,*)& "All models have irreparable overlaps. Trying randoms starts." iranconf=1 i_model=nmodel_start+1 goto 122 endif else ! Remove SC overlaps if requested if (overlapsc) then write (iout,*) 'Calling OVERLAP_SC' call overlap_sc(fail) if (fail) then write (iout,*)& "Failed to remove overlap" endif endif if (me.eq.king.or..not.out1file) then write (iout,*) "Energies after removing overlaps" call etotal(energia(0)) call enerprint(energia(0)) endif endif ! 8/22/17 AL Minimize initial structure if (dccart) then if (me.eq.king.or..not.out1file) write(iout,*)& 'Minimizing initial PDB structure: Calling MINIM_DC' call minim_dc(etot,iretcode,nfun) #ifdef LBFGS if (me.eq.king.or..not.out1file)& write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun #endif else call geom_to_var(nvar,varia) if(me.eq.king.or..not.out1file) write (iout,*)& 'Minimizing initial PDB structure: Calling MINIMIZE.' call minimize(etot,varia,iretcode,nfun) call var_to_geom(nvar,varia) #ifdef LBFGS if (me.eq.king.or..not.out1file)& write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun if(me.eq.king.or..not.out1file)& write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun #else if (me.eq.king.or..not.out1file)& write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun if(me.eq.king.or..not.out1file)& write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun #endif endif endif if (nmodel_start.gt.0 .and. me.eq.king) then write (iout,'(a)') "Task Starting model" do i=0,nodes-1 if (i_start_models(i).gt.nmodel_start) then write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE" else write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i)) & (:ilen(pdbfiles_chomo(i_start_models(i)))) endif enddo endif endif call chainbuild_cart call kinetic(EK) write(iout,*) "just after kinetic" call cartprint if (tbf) then call verlet_bath endif kinetic_T=2.0d0/(dimen3*Rb)*EK if(me.eq.king.or..not.out1file)then write(iout,*) "just after verlet_bath" call cartprint call intout endif #ifdef MPI tt0=MPI_Wtime() #else tt0=tcpu() #endif call zerograd write(iout,*) "before ETOTAL" call etotal(potEcomp) if (large) call enerprint(potEcomp) #ifdef TIMING_ENE #ifdef MPI t_etotal=t_etotal+MPI_Wtime()-tt0 #else t_etotal=t_etotal+tcpu()-tt0 #endif #endif potE=potEcomp(0) call cartgrad write(iout,*) "before lagrangian" call lagrangian write(iout,*) "before max_accel" call max_accel if (amax*d_time .gt. dvmax) then d_time=d_time*dvmax/amax if(me.eq.king.or..not.out1file) write (iout,*) & "Time step reduced to",d_time,& " because of too large initial acceleration." endif if(me.eq.king.or..not.out1file)then write(iout,*) "Potential energy and its components" call enerprint(potEcomp) ! write(iout,*) (potEcomp(i),i=0,n_ene) endif potE=potEcomp(0)-potEcomp(51) totE=EK+potE itime=0 itime_mat=itime if (ntwe.ne.0) call statout(itime) if(me.eq.king.or..not.out1file) & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", & " Kinetic energy",EK," Potential energy",potE, & " Total energy",totE," Maximum acceleration ", & amax if (large) then write (iout,*) "Initial coordinates" do i=1,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),& (c(j,i+nres),j=1,3) enddo write (iout,*) "Initial dC" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),& (dc(j,i+nres),j=1,3) enddo write (iout,*) "Initial velocities" write (iout,"(13x,' backbone ',23x,' side chain')") do i=0,nres write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo write (iout,*) "Initial accelerations" do i=0,nres ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),& (d_a(j,i+nres),j=1,3) enddo endif do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3) enddo if (RESPA) then #ifdef MPI tt0 =MPI_Wtime() #else tt0 = tcpu() #endif call zerograd call etotal_short(energia_short) if (large) call enerprint(potEcomp) #ifdef TIMING_ENE #ifdef MPI t_eshort=t_eshort+MPI_Wtime()-tt0 #else t_eshort=t_eshort+tcpu()-tt0 #endif #endif call cartgrad call lagrangian if(.not.out1file .and. large) then write (iout,*) "energia_long",energia_long(0),& " energia_short",energia_short(0),& " total",energia_long(0)+energia_short(0) write (iout,*) "Initial fast-force accelerations" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),& (d_a(j,i+nres),j=1,3) enddo endif ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array do i=0,2*nres do j=1,3 d_a_short(j,i)=d_a(j,i) enddo enddo #ifdef MPI tt0=MPI_Wtime() #else tt0=tcpu() #endif call zerograd call etotal_long(energia_long) if (large) call enerprint(potEcomp) #ifdef TIMING_ENE #ifdef MPI t_elong=t_elong+MPI_Wtime()-tt0 #else t_elong=t_elong+tcpu()-tt0 #endif #endif call cartgrad call lagrangian if(.not.out1file .and. large) then write (iout,*) "energia_long",energia_long(0) write (iout,*) "Initial slow-force accelerations" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),& (d_a(j,i+nres),j=1,3) enddo endif #ifdef MPI t_enegrad=t_enegrad+MPI_Wtime()-tt0 #else t_enegrad=t_enegrad+tcpu()-tt0 #endif endif return end subroutine init_MD !----------------------------------------------------------------------------- subroutine random_vel ! implicit real*8 (a-h,o-z) use energy_data use random, only:anorm_distr use MD_data ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' real(kind=8) :: xv,sigv,lowb,highb ,Ek1 #ifdef FIVEDIAG integer ichain,n,innt,inct,ibeg,ierr,innt_org real(kind=8) ,allocatable, dimension(:):: work integer,allocatable,dimension(:) :: iwork ! double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),& ! Gvec(maxres2_chain,maxres2_chain) ! common /przechowalnia/Ghalf,Geigen,Gvec !#ifdef DEBUG ! double precision inertia(maxres2_chain,maxres2_chain) !#endif #endif !#define DEBUG #ifdef FIVEDIAG real(kind=8) ,allocatable, dimension(:) :: xsolv,DML,rs real(kind=8) :: sumx,Ek2,Ek3,aux,masinv #ifdef DEBUG real(kind=8) ,allocatable, dimension(:) :: rsold real (kind=8),allocatable,dimension(:,:) :: matold,inertia integer :: iti #endif #endif integer :: i,j,ii,k,mark,imark,mnum,nres2 integer(kind=8) :: ind ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m !#undef DEBUG ! First generate velocities in the eigenspace of the G matrix ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3 ! call flush(iout) #ifdef FIVEDIAG if(.not.allocated(work)) then allocate(work(48*nres)) allocate(iwork(6*nres)) endif print *,"IN RANDOM VEL" nres2=2*nres ! print *,size(ghalf) #undef DEBUG #ifdef DEBUG write (iout,*) "Random_vel, fivediag" flush(iout) allocate(inertia(2*nres,2*nres)) #endif d_t=0.0d0 Ek2=0.0d0 EK=0.0d0 Ek3=0.0d0 #ifdef DEBUG write(iout,*), "nchain",nchain #endif do ichain=1,nchain ind=0 ! if(.not.allocated(ghalf)) print *,"COCO" ! if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) ! ghalf=0.0d0 n=dimen_chain(ichain) innt=iposd_chain(ichain) ! innt_org= innt_org=chain_border(1,ichain) if ((molnum(innt_org).eq.5).or.(molnum(innt_org).eq.4)) go to 137 if(.not.allocated(ghalf)) print *,"COCO" if(.not.allocated(Ghalf)) allocate(Ghalf(1300*(1300+1)/2)) ghalf=0.0d0 inct=innt+n-1 #ifdef DEBUG write (iout,*) "Chain",ichain," n",n," start",innt do i=innt,inct if (i.lt.inct-1) then write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),& DU2orig(i) else if (i.eq.inct-1) then write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i) else write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i) endif enddo #endif ghalf(ind+1)=dmorig(innt) ghalf(ind+2)=du1orig(innt) ghalf(ind+3)=dmorig(innt+1) ind=ind+3 do i=3,n ind=ind+i-3 write (iout,*) "i",i," ind",ind," indu2",innt+i-2,& " indu1",innt+i-1," indm",innt+i ghalf(ind+1)=du2orig(innt-1+i-2) ghalf(ind+2)=du1orig(innt-1+i-1) ghalf(ind+3)=dmorig(innt-1+i) !c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2, !c & "DU1",innt-1+i-1,"DM ",innt-1+i ind=ind+3 enddo #ifdef DEBUG ind=0 do i=1,n do j=1,i ind=ind+1 inertia(i,j)=ghalf(ind) inertia(j,i)=ghalf(ind) enddo enddo #endif #ifdef DEBUG write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2 write (iout,*) "Five-diagonal inertia matrix, lower triangle" ! call matoutr(n,ghalf) #endif call gldiag(nres*2,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork) if (large) then write (iout,'(//a,i3)')& "Eigenvectors and eigenvalues of the G matrix chain",ichain call eigout(n,n,nres*2,nres*2,Gvec,Geigen) endif #ifdef DIAGCHECK !c check diagonalization do i=1,n do j=1,n aux=0.0d0 do k=1,n do l=1,n aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l) enddo enddo if (i.eq.j) then write (iout,*) i,j,aux,geigen(i) else write (iout,*) i,j,aux endif enddo enddo #endif 137 continue write(iout,*) "HERE,",n,innt innt_org=chain_border(1,ichain) xv=0.0d0 ii=0 do i=1,n do k=1,3 ii=ii+1 mnum=molnum(innt_org) if (molnum(innt_org).ge.4) geigen(i)=3.0/msc(itype(innt_org+i-1,mnum),mnum) ! if (molnum(innt).eq.5) write(iout,*) "typ",i,innt-1+i,itype(innt+i-1,5) sigv=dsqrt((Rb*t_bath)/geigen(i)) lowb=-5*sigv highb=5*sigv d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb) EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2 write (iout,*) "i",i," ii",ii," geigen",geigen(i), & " d_t_work_new",d_t_work_new(ii),innt_org+i-1 enddo enddo if (molnum(innt_org).ge.4) then mnum=molnum(innt_org) do k=1,3 do i=1,n ind=(i-1)*3+k d_t_work(ind)=0.0d0 masinv=1.0d0/msc(itype(innt_org+i-1,mnum),mnum) d_t_work(ind)=d_t_work(ind)& +masinv*d_t_work_new((i-1)*3+k) enddo enddo else do k=1,3 do i=1,n ind=(i-1)*3+k d_t_work(ind)=0.0d0 do j=1,n d_t_work(ind)=d_t_work(ind)& +Gvec(i,j)*d_t_work_new((j-1)*3+k) enddo enddo enddo endif #ifdef DEBUG aux=0.0d0 do k=1,3 do i=1,n do j=1,n aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k) enddo enddo enddo Ek3=Ek3+aux/2 #endif !c Transfer to the d_t vector innt=chain_border(1,ichain) inct=chain_border(2,ichain) ind=0 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct do i=innt,inct do j=1,3 ind=ind+1 d_t(j,i)=d_t_work(ind) enddo mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.le.2) then do j=1,3 ind=ind+1 d_t(j,i+nres)=d_t_work(ind) enddo endif enddo enddo if (large) then write (iout,*) write (iout,*) "Random velocities in the Calpha,SC space" do i=1,nres mnum=molnum(i) write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')& restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3) enddo endif call kinetic_CASC(Ek1) ! ! Transform the velocities to virtual-bond space ! #define WLOS #ifdef WLOS if (nnt.eq.1) then d_t(:,0)=d_t(:,1) endif do i=1,nres mnum=molnum(i) if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or.mnum.ge.3) then do j=1,3 d_t(j,i)=d_t(j,i+1)-d_t(j,i) enddo else do j=1,3 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i) d_t(j,i)=d_t(j,i+1)-d_t(j,i) enddo end if enddo d_t(:,nres)=0.0d0 d_t(:,nct)=0.0d0 d_t(:,2*nres)=0.0d0 if (nnt.gt.1) then d_t(:,0)=d_t(:,1) d_t(:,1)=0.0d0 endif if (large) then write (iout,*) write (iout,*) "Random vel after 1st transf the Calpha,SC space" write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3) do i=1,nres mnum=molnum(i) write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')& restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3) enddo endif !c d_a(:,0)=d_a(:,1) !c d_a(:,1)=0.0d0 !c write (iout,*) "Shifting accelerations" do ichain=2,nchain write(iout,*) "nchain",ichain,chain_border1(1,ichain),molnum(chain_border1(1,ichain)) if (molnum(chain_border1(1,ichain)+1).eq.5) cycle !c write (iout,*) "ichain",chain_border1(1,ichain)-1, !c & chain_border1(1,ichain) d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain)) d_t(:,chain_border1(1,ichain))=0.0d0 enddo !c write (iout,*) "Adding accelerations" do ichain=2,nchain if (molnum(chain_border1(1,ichain)+1).eq.5) cycle !c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1, !c & chain_border(2,ichain-1) d_t(:,chain_border1(1,ichain)-1)=& d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1)) d_t(:,chain_border(2,ichain-1))=0.0d0 enddo do ichain=2,nchain write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,& chain_border(2,ichain-1) if (molnum(chain_border1(1,ichain)+1).eq.5) cycle d_t(:,chain_border1(1,ichain)-1)=& d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1)) d_t(:,chain_border(2,ichain-1))=0.0d0 enddo if (large) then write (iout,*) write (iout,*) "Random vel after 2nd transf the Calpha,SC space" write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3) do i=1,nres mnum=molnum(i) write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')& restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3) enddo endif #else ibeg=0 !c do j=1,3 !c d_t(j,0)=d_t(j,nnt) !c enddo do ichain=1,nchain innt=chain_border(1,ichain) inct=chain_border(2,ichain) !c write (iout,*) "ichain",ichain," innt",innt," inct",inct !c write (iout,*) "ibeg",ibeg do j=1,3 d_t(j,ibeg)=d_t(j,innt) enddo ibeg=inct+1 do i=innt,inct mnum=molnum(i) if (iabs(itype(i,1).eq.10).or.mnum.ge.3) then !c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3) do j=1,3 d_t(j,i)=d_t(j,i+1)-d_t(j,i) enddo else do j=1,3 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i) d_t(j,i)=d_t(j,i+1)-d_t(j,i) enddo end if enddo enddo #endif if (large) then write (iout,*) write (iout,*)& "Random velocities in the virtual-bond-vector space" write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3) do i=1,nres write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')& restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3) enddo write (iout,*) write (iout,*) "Kinetic energy from inertia matrix eigenvalues",& Ek write (iout,*)& "Kinetic temperatures from inertia matrix eigenvalues",& 2*Ek/(3*dimen*Rb) #ifdef DEBUG write (iout,*) "Kinetic energy from inertia matrix",Ek3 write (iout,*) "Kinetic temperatures from inertia",& 2*Ek3/(3*dimen*Rb) #endif write (iout,*) "Kinetic energy from velocities in CA-SC space",& Ek1 write (iout,*)& "Kinetic temperatures from velovities in CA-SC space",& 2*Ek1/(3*dimen*Rb) call kinetic(Ek1) write (iout,*)& "Kinetic energy from virtual-bond-vector velocities",Ek1 write (iout,*)& "Kinetic temperature from virtual-bond-vector velocities ",& 2*Ek1/(dimen3*Rb) endif #else xv=0.0d0 ii=0 do i=1,dimen do k=1,3 ii=ii+1 sigv=dsqrt((Rb*t_bath)/geigen(i)) lowb=-5*sigv highb=5*sigv d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb) #ifdef DEBUG write (iout,*) "i",i," ii",ii," geigen",geigen(i),& " d_t_work_new",d_t_work_new(ii) #endif enddo enddo #ifdef DEBUG ! diagnostics Ek1=0.0d0 ii=0 do i=1,dimen do k=1,3 ii=ii+1 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2 enddo enddo write (iout,*) "Ek from eigenvectors",Ek1 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) ! end diagnostics #endif do k=0,2 do i=1,dimen ind=(i-1)*3+k+1 d_t_work(ind)=0.0d0 do j=1,dimen d_t_work(ind)=d_t_work(ind) & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1) enddo ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind) ! call flush(iout) enddo enddo ! Transfer to the d_t vector do j=1,3 d_t(j,0)=d_t_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 ind=ind+1 d_t(j,i)=d_t_work(ind) enddo enddo do i=nnt,nct mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then do j=1,3 ind=ind+1 d_t(j,i+nres)=d_t_work(ind) enddo endif enddo #endif ! call kinetic(EK) ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",& ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1 ! call flush(iout) ! write(iout,*) "end init MD" #undef DEBUG return end subroutine random_vel !----------------------------------------------------------------------------- #ifndef LANG0 subroutine sd_verlet_p_setup ! Sets up the parameters of stochastic Verlet algorithm ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' use control, only: tcpu use control_data #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres real(kind=8) :: pterm,vterm,rho,rhoc,vsig real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,& prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres logical :: lprn = .false. real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,& gdt9,psig,tt0 integer :: i,maxres2 #ifdef MPI tt0 = MPI_Wtime() #else tt0 = tcpu() #endif ! ! AL 8/17/04 Code adapted from tinker ! ! Get the frictional and random terms for stochastic dynamics in the ! eigenspace of mass-scaled UNRES friction matrix ! maxres2=2*nres do i = 1, dimen gdt = fricgam(i) * d_time ! ! Stochastic dynamics reduces to simple MD for zero friction ! if (gdt .le. zero) then pfric_vec(i) = 1.0d0 vfric_vec(i) = d_time afric_vec(i) = 0.5d0 * d_time * d_time prand_vec(i) = 0.0d0 vrand_vec1(i) = 0.0d0 vrand_vec2(i) = 0.0d0 ! ! Analytical expressions when friction coefficient is large ! else if (gdt .ge. gdt_radius) then egdt = dexp(-gdt) pfric_vec(i) = egdt vfric_vec(i) = (1.0d0-egdt) / fricgam(i) afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i) pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt vterm = 1.0d0 - egdt**2 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm) ! ! Use series expansions when friction coefficient is small ! else gdt2 = gdt * gdt gdt3 = gdt * gdt2 gdt4 = gdt2 * gdt2 gdt5 = gdt2 * gdt3 gdt6 = gdt3 * gdt3 gdt7 = gdt3 * gdt4 gdt8 = gdt4 * gdt4 gdt9 = gdt4 * gdt5 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 & - gdt5/120.0d0 + gdt6/720.0d0 & - gdt7/5040.0d0 + gdt8/40320.0d0 & - gdt9/362880.0d0) / fricgam(i)**2 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i) pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i) pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 & + 127.0d0*gdt9/90720.0d0 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 & - 17.0d0*gdt2/1280.0d0 & + 17.0d0*gdt3/6144.0d0 & + 40967.0d0*gdt4/34406400.0d0 & - 57203.0d0*gdt5/275251200.0d0 & - 1429487.0d0*gdt6/13212057600.0d0) end if ! ! Compute the scaling factors of random terms for the nonzero friction case ! ktm = 0.5d0*d_time/fricgam(i) psig = dsqrt(ktm*pterm) / fricgam(i) vsig = dsqrt(ktm*vterm) rhoc = dsqrt(1.0d0 - rho*rho) prand_vec(i) = psig vrand_vec1(i) = vsig * rho vrand_vec2(i) = vsig * rhoc end if end do if (lprn) then write (iout,*) & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",& " vrand_vec2" do i=1,dimen write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),& afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i) enddo endif ! ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables ! #ifndef LANG0 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat) call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat) call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat) call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat) call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1) call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2) #endif #ifdef MPI t_sdsetup=t_sdsetup+MPI_Wtime() #else t_sdsetup=t_sdsetup+tcpu()-tt0 #endif return end subroutine sd_verlet_p_setup !----------------------------------------------------------------------------- subroutine eigtransf1(n,ndim,ab,d,c) !el implicit none integer :: n,ndim real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim) integer :: i,j,k do i=1,n do j=1,n c(i,j)=0.0d0 do k=1,n c(i,j)=c(i,j)+ab(k,j,i)*d(k) enddo enddo enddo return end subroutine eigtransf1 !----------------------------------------------------------------------------- subroutine eigtransf(n,ndim,a,b,d,c) !el implicit none integer :: n,ndim real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim) integer :: i,j,k do i=1,n do j=1,n c(i,j)=0.0d0 do k=1,n c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k) enddo enddo enddo return end subroutine eigtransf !----------------------------------------------------------------------------- subroutine sd_verlet1 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec logical :: lprn = .false. real(kind=8) :: ddt1,ddt2 integer :: i,j,ind,inres ! write (iout,*) "dc_old" ! do i=0,nres ! write (iout,'(i5,3f10.5,5x,3f10.5)') ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3) ! enddo do j=1,3 dc_work(j)=dc_old(j,0) d_t_work(j)=d_t_old(j,0) d_a_work(j)=d_a_old(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 dc_work(ind+j)=dc_old(j,i) d_t_work(ind+j)=d_t_old(j,i) d_a_work(ind+j)=d_a_old(j,i) enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) d_a_work(ind+j)=d_a_old(j,i+nres) enddo ind=ind+3 endif enddo #ifndef LANG0 if (lprn) then write (iout,*) & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",& " vrand_mat2" do i=1,dimen do j=1,dimen write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),& vfric_mat(i,j),afric_mat(i,j),& prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j) enddo enddo endif do i=1,dimen ddt1=0.0d0 ddt2=0.0d0 do j=1,dimen dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j) ddt1=ddt1+pfric_mat(i,j)*d_t_work(j) ddt2=ddt2+vfric_mat(i,j)*d_a_work(j) enddo d_t_work_new(i)=ddt1+0.5d0*ddt2 d_t_work(i)=ddt1+ddt2 enddo #endif do j=1,3 dc(j,0)=dc_work(j) d_t(j,0)=d_t_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 dc(j,i)=dc_work(ind+j) d_t(j,i)=d_t_work(ind+j) enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then inres=i+nres do j=1,3 dc(j,inres)=dc_work(ind+j) d_t(j,inres)=d_t_work(ind+j) enddo ind=ind+3 endif enddo return end subroutine sd_verlet1 !----------------------------------------------------------------------------- subroutine sd_verlet2 ! Calculating the adjusted velocities for accelerations use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec ! real(kind=8) :: ddt1,ddt2 integer :: i,j,ind,inres ! Compute the stochastic forces which contribute to velocity change ! call stochastic_force(stochforcvecV) #ifndef LANG0 do i=1,dimen ddt1=0.0d0 ddt2=0.0d0 do j=1,dimen ddt1=ddt1+vfric_mat(i,j)*d_a_work(j) ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ & vrand_mat2(i,j)*stochforcvecV(j) enddo d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2 enddo #endif do j=1,3 d_t(j,0)=d_t_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_work(ind+j) enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j) enddo ind=ind+3 endif enddo return end subroutine sd_verlet2 !----------------------------------------------------------------------------- subroutine sd_verlet_ciccotti_setup ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's ! version ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' use control, only: tcpu use control_data #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres real(kind=8) :: pterm,vterm,rho,rhoc,vsig real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,& prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres logical :: lprn = .false. real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0 real(kind=8) :: ktm,gdt,egdt,tt0 integer :: i,maxres2 #ifdef MPI tt0 = MPI_Wtime() #else tt0 = tcpu() #endif ! ! AL 8/17/04 Code adapted from tinker ! ! Get the frictional and random terms for stochastic dynamics in the ! eigenspace of mass-scaled UNRES friction matrix ! maxres2=2*nres do i = 1, dimen write (iout,*) "i",i," fricgam",fricgam(i) gdt = fricgam(i) * d_time ! ! Stochastic dynamics reduces to simple MD for zero friction ! if (gdt .le. zero) then pfric_vec(i) = 1.0d0 vfric_vec(i) = d_time afric_vec(i) = 0.5d0*d_time*d_time prand_vec(i) = afric_vec(i) vrand_vec2(i) = vfric_vec(i) ! ! Analytical expressions when friction coefficient is large ! else egdt = dexp(-gdt) pfric_vec(i) = egdt vfric_vec(i) = dexp(-0.5d0*gdt)*d_time afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time prand_vec(i) = afric_vec(i) vrand_vec2(i) = vfric_vec(i) ! ! Compute the scaling factors of random terms for the nonzero friction case ! ! ktm = 0.5d0*d_time/fricgam(i) ! psig = dsqrt(ktm*pterm) / fricgam(i) ! vsig = dsqrt(ktm*vterm) ! prand_vec(i) = psig*afric_vec(i) ! vrand_vec2(i) = vsig*vfric_vec(i) end if end do if (lprn) then write (iout,*) & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",& " vrand_vec2" do i=1,dimen write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),& afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i) enddo endif ! ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables ! call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat) call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat) call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat) call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat) call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2) #ifdef MPI t_sdsetup=t_sdsetup+MPI_Wtime() #else t_sdsetup=t_sdsetup+tcpu()-tt0 #endif return end subroutine sd_verlet_ciccotti_setup !----------------------------------------------------------------------------- subroutine sd_verlet1_ciccotti ! Applying stochastic velocity Verlet algorithm - step 1 to velocities ! implicit real*8 (a-h,o-z) use energy_data ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec logical :: lprn = .false. real(kind=8) :: ddt1,ddt2 integer :: i,j,ind,inres ! write (iout,*) "dc_old" ! do i=0,nres ! write (iout,'(i5,3f10.5,5x,3f10.5)') ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3) ! enddo do j=1,3 dc_work(j)=dc_old(j,0) d_t_work(j)=d_t_old(j,0) d_a_work(j)=d_a_old(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 dc_work(ind+j)=dc_old(j,i) d_t_work(ind+j)=d_t_old(j,i) d_a_work(ind+j)=d_a_old(j,i) enddo ind=ind+3 enddo do i=nnt,nct if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) d_a_work(ind+j)=d_a_old(j,i+nres) enddo ind=ind+3 endif enddo #ifndef LANG0 if (lprn) then write (iout,*) & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",& " vrand_mat2" do i=1,dimen do j=1,dimen write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),& vfric_mat(i,j),afric_mat(i,j),& prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j) enddo enddo endif do i=1,dimen ddt1=0.0d0 ddt2=0.0d0 do j=1,dimen dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j) ddt1=ddt1+pfric_mat(i,j)*d_t_work(j) ddt2=ddt2+vfric_mat(i,j)*d_a_work(j) enddo d_t_work_new(i)=ddt1+0.5d0*ddt2 d_t_work(i)=ddt1+ddt2 enddo #endif do j=1,3 dc(j,0)=dc_work(j) d_t(j,0)=d_t_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 dc(j,i)=dc_work(ind+j) d_t(j,i)=d_t_work(ind+j) enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then inres=i+nres do j=1,3 dc(j,inres)=dc_work(ind+j) d_t(j,inres)=d_t_work(ind+j) enddo ind=ind+3 endif enddo return end subroutine sd_verlet1_ciccotti !----------------------------------------------------------------------------- subroutine sd_verlet2_ciccotti ! Calculating the adjusted velocities for accelerations use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec real(kind=8) :: ddt1,ddt2 integer :: i,j,ind,inres ! ! Compute the stochastic forces which contribute to velocity change ! call stochastic_force(stochforcvecV) #ifndef LANG0 do i=1,dimen ddt1=0.0d0 ddt2=0.0d0 do j=1,dimen ddt1=ddt1+vfric_mat(i,j)*d_a_work(j) ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j) ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j) enddo d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2 enddo #endif do j=1,3 d_t(j,0)=d_t_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 d_t(j,i)=d_t_work(ind+j) enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j) enddo ind=ind+3 endif enddo return end subroutine sd_verlet2_ciccotti #endif !----------------------------------------------------------------------------- ! moments.f !----------------------------------------------------------------------------- #ifdef FIVEDIAG subroutine inertia_tensor use comm_gucio use energy_data real(kind=8) Im(3,3),Imcp(3,3),pr(3),M_SC,& eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),& vpp(3,0:MAXRES),vs_p(3),pr1(3,3),& pr2(3,3),pp(3),incr(3),v(3),mag,mag2,M_PEP integer iti,inres,i,j,k,mnum,mnum1 do i=1,3 do j=1,3 Im(i,j)=0.0d0 pr1(i,j)=0.0d0 pr2(i,j)=0.0d0 enddo L(i)=0.0d0 cm(i)=0.0d0 vrot(i)=0.0d0 enddo M_PEP=0.0d0 !c caulating the center of the mass of the protein do i=nnt,nct-1 mnum=molnum(i) mnum1=molnum(i+1) if (itype(i,mnum).eq.ntyp1_molec(mnum)& .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) if (mnum.ge.5) mp(mnum)=0.0d0 M_PEP=M_PEP+mp(mnum) do j=1,3 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum) enddo enddo ! do j=1,3 ! cm(j)=mp*cm(j) ! enddo M_SC=0.0d0 do i=nnt,nct mnum=molnum(i) mnum1=molnum(i+1) iti=iabs(itype(i,mnum)) if (iti.eq.ntyp1_molec(mnum)) cycle M_SC=M_SC+msc(iabs(iti),mnum) inres=i+nres do j=1,3 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres) enddo enddo do j=1,3 cm(j)=cm(j)/(M_SC+M_PEP) enddo do i=nnt,nct-1 mnum=molnum(i) mnum1=molnum(i+1) ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) if (mnum.ge.5) mp(mnum)=0.0d0 if (itype(i,mnum).eq.ntyp1_molec(mnum)& .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3)) Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2) Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3) Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3) Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1)) Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2)) enddo do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) if (iti.eq.ntyp1_molec(mnum)) cycle inres=i+nres do j=1,3 pr(j)=c(j,inres)-cm(j) enddo Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3)) Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2) Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3) Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3) Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1)) Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2)) enddo do i=nnt,nct-1 mnum=molnum(i) mnum1=molnum(i+1) if (itype(i,mnum).eq.ntyp1_molec(mnum)& .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))*& vbld(i+1)*vbld(i+1)*0.25d0 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))*& vbld(i+1)*vbld(i+1)*0.25d0 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))*& vbld(i+1)*vbld(i+1)*0.25d0 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))*& vbld(i+1)*vbld(i+1)*0.25d0 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))*& vbld(i+1)*vbld(i+1)*0.25d0 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*& vbld(i+1)*vbld(i+1)*0.25d0 enddo do i=nnt,nct mnum=molnum(i) mnum1=molnum(i+1) iti=iabs(itype(i,mnum)) if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum).and.mnum.le.2) then inres=i+nres Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)*& dc_norm(1,inres))*vbld(inres)*vbld(inres) Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)*& dc_norm(2,inres))*vbld(inres)*vbld(inres) Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)*& dc_norm(3,inres))*vbld(inres)*vbld(inres) Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)*& dc_norm(3,inres))*vbld(inres)*vbld(inres) Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)*& dc_norm(2,inres))*vbld(inres)*vbld(inres) Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*& dc_norm(3,inres))*vbld(inres)*vbld(inres) endif enddo call angmom(cm,L) Im(2,1)=Im(1,2) Im(3,1)=Im(1,3) Im(3,2)=Im(2,3) !c Copng the Im matrix for the djacob subroutine do i=1,3 do j=1,3 Imcp(i,j)=Im(i,j) Id(i,j)=0.0d0 enddo enddo !c Finding the eigenvectors and eignvalues of the inertia tensor call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval) do i=1,3 if (dabs(eigval(i)).gt.1.0d-15) then Id(i,i)=1.0d0/eigval(i) else Id(i,i)=0.0d0 endif enddo do i=1,3 do j=1,3 Imcp(i,j)=eigvec(j,i) enddo enddo do i=1,3 do j=1,3 do k=1,3 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j) enddo enddo enddo do i=1,3 do j=1,3 do k=1,3 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j) enddo enddo enddo !c Calculating the total rotational velocity of the molecule do i=1,3 do j=1,3 vrot(i)=vrot(i)+pr2(i,j)*L(j) enddo enddo do i=nnt,nct-1 mnum=molnum(i) mnum1=molnum(i+1) ! write (iout,*) itype(i+1,mnum1),itype(i,mnum) if (itype(i+1,mnum1).ne.ntyp1_molec(mnum1) & .and. itype(i,mnum).eq.ntyp1_molec(mnum) .or.& itype(i,mnum).ne.ntyp1_molec(mnum) & .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle call vecpr(vrot(1),dc(1,i),vp) do j=1,3 d_t(j,i)=d_t(j,i)-vp(j) enddo enddo do i=nnt,nct mnum=molnum(i) if(itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)& .and.mnum.le.2) then inres=i+nres call vecpr(vrot(1),dc(1,inres),vp) do j=1,3 d_t(j,inres)=d_t(j,inres)-vp(j) enddo endif enddo call angmom(cm,L) return end subroutine inertia_tensor !------------------------------------------------------------ subroutine angmom(cm,L) implicit none double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),& pp(3),mscab integer iti,inres,i,j,mnum,mnum1 !c Calculate the angular momentum do j=1,3 L(j)=0.0d0 enddo do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct-1 mnum=molnum(i) mnum1=molnum(i+1) ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) if (mnum.ge.5) mp(mnum)=0.0d0 if (itype(i,mnum).eq.ntyp1_molec(mnum)& .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) enddo do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo call vecpr(pr(1),v(1),vp) do j=1,3 L(j)=L(j)+mp(mnum)*vp(j) enddo do j=1,3 pr(j)=0.5d0*dc(j,i) pp(j)=0.5d0*d_t(j,i) enddo call vecpr(pr(1),pp(1),vp) do j=1,3 L(j)=L(j)+Ip(mnum)*vp(j) enddo enddo do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) inres=i+nres do j=1,3 pr(j)=c(j,inres)-cm(j) enddo if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)& .and.mnum.le.2) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo else do j=1,3 v(j)=incr(j) enddo endif call vecpr(pr(1),v(1),vp) !c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3), !c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3) ! if (mnum.gt.4) then ! mscab=0.0d0 ! else mscab=msc(iti,mnum) ! endif do j=1,3 L(j)=L(j)+mscab*vp(j) enddo !c write (iout,*) "L",(l(j),j=1,3) if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)& .and.mnum.le.2) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo call vecpr(dc(1,inres),d_t(1,inres),vp) do j=1,3 L(j)=L(j)+Isc(iti,mnum)*vp(j) enddo endif do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo enddo return end subroutine angmom !--------------------------------------------------------------- subroutine vcm_vel(vcm) double precision vcm(3),vv(3),summas,amas integer i,j,mnum,mnum1 do j=1,3 vcm(j)=0.0d0 vv(j)=d_t(j,0) enddo summas=0.0d0 do i=nnt,nct mnum=molnum(i) if ((mnum.gt.5).or.(mnum.eq.3))& mp(mnum)=msc(itype(i,mnum),mnum) if (i.lt.nct) then summas=summas+mp(mnum) do j=1,3 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i)) enddo endif ! if (mnum.ne.4) then amas=msc(iabs(itype(i,mnum)),mnum) ! else ! amas=0.0d0 ! endif ! amas=msc(iabs(itype(i))) summas=summas+amas ! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then do j=1,3 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres)) enddo else do j=1,3 vcm(j)=vcm(j)+amas*vv(j) enddo endif do j=1,3 vv(j)=vv(j)+d_t(j,i) enddo enddo write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas do j=1,3 vcm(j)=vcm(j)/summas enddo return end subroutine vcm_vel #else subroutine inertia_tensor ! Calculating the intertia tensor for the entire protein in order to ! remove the perpendicular components of velocity matrix which cause ! the molecule to rotate. use comm_gucio use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot real(kind=8) :: M_SC,mag,mag2,M_PEP real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES) real(kind=8),dimension(3) :: vs_p,pp,incr,v real(kind=8),dimension(3,3) :: pr1,pr2 !el common /gucio/ cm integer :: iti,inres,i,j,k,mnum do i=1,3 do j=1,3 Im(i,j)=0.0d0 pr1(i,j)=0.0d0 pr2(i,j)=0.0d0 enddo L(i)=0.0d0 cm(i)=0.0d0 vrot(i)=0.0d0 enddo ! calculating the center of the mass of the protein M_PEP=0.0d0 do i=nnt,nct-1 mnum=molnum(i) ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) write(iout,*) "WTF",itype(i,mnum),i,mnum,mp(mnum) ! if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle M_PEP=M_PEP+mp(mnum) do j=1,3 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum) enddo enddo ! do j=1,3 ! cm(j)=mp(1)*cm(j) ! enddo M_SC=0.0d0 do i=nnt,nct mnum=molnum(i) ! if (mnum.ge.5) cycle iti=iabs(itype(i,mnum)) M_SC=M_SC+msc(iabs(iti),mnum) inres=i+nres if (mnum.ge.4) inres=i do j=1,3 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres) enddo enddo do j=1,3 cm(j)=cm(j)/(M_SC+M_PEP) enddo ! write(iout,*) "Center of mass:",cm do i=nnt,nct-1 mnum=molnum(i) ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3)) Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2) Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3) Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3) Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1)) Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2)) enddo ! write(iout,*) "The angular momentum before msc add" ! do i=1,3 ! write (iout,*) (Im(i,j),j=1,3) ! enddo do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) ! if (mnum.ge.5) cycle inres=i+nres if (mnum.ge.4) inres=i do j=1,3 pr(j)=c(j,inres)-cm(j) enddo Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3)) Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2) Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3) Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3) Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1)) Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2)) enddo ! write(iout,*) "The angular momentum before Ip add" ! do i=1,3 ! write (iout,*) (Im(i,j),j=1,3) ! enddo do i=nnt,nct-1 mnum=molnum(i) Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 enddo ! write(iout,*) "The angular momentum before Isc add" ! do i=1,3 ! write (iout,*) (Im(i,j),j=1,3) ! enddo do i=nnt,nct mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then iti=iabs(itype(i,mnum)) inres=i+nres Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* & dc_norm(1,inres))*vbld(inres)*vbld(inres) Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* & dc_norm(2,inres))*vbld(inres)*vbld(inres) Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* & dc_norm(2,inres))*vbld(inres)*vbld(inres) Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) endif enddo ! write(iout,*) "The angular momentum before agnom:" ! do i=1,3 ! write (iout,*) (Im(i,j),j=1,3) ! enddo call angmom(cm,L) ! write(iout,*) "The angular momentum before adjustment:" ! write(iout,*) (L(j),j=1,3) ! do i=1,3 ! write (iout,*) (Im(i,j),j=1,3) ! enddo Im(2,1)=Im(1,2) Im(3,1)=Im(1,3) Im(3,2)=Im(2,3) ! Copying the Im matrix for the djacob subroutine do i=1,3 do j=1,3 Imcp(i,j)=Im(i,j) Id(i,j)=0.0d0 enddo enddo ! Finding the eigenvectors and eignvalues of the inertia tensor call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval) ! write (iout,*) "Eigenvalues & Eigenvectors" ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3) ! write (iout,*) ! do i=1,3 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3) ! enddo ! Constructing the diagonalized matrix do i=1,3 if (dabs(eigval(i)).gt.1.0d-15) then Id(i,i)=1.0d0/eigval(i) else Id(i,i)=0.0d0 endif enddo do i=1,3 do j=1,3 Imcp(i,j)=eigvec(j,i) enddo enddo do i=1,3 do j=1,3 do k=1,3 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j) enddo enddo enddo do i=1,3 do j=1,3 do k=1,3 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j) enddo enddo enddo ! Calculating the total rotational velocity of the molecule do i=1,3 do j=1,3 vrot(i)=vrot(i)+pr2(i,j)*L(j) enddo enddo ! Resetting the velocities do i=nnt,nct-1 call vecpr(vrot(1),dc(1,i),vp) do j=1,3 ! print *,"HERE2",d_t(j,i),vp(j) d_t(j,i)=d_t(j,i)-vp(j) ! print *,"HERE2",d_t(j,i),vp(j) enddo enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then inres=i+nres call vecpr(vrot(1),dc(1,inres),vp) do j=1,3 d_t(j,inres)=d_t(j,inres)-vp(j) enddo endif enddo call angmom(cm,L) ! write(iout,*) "The angular momentum after adjustment:" ! write(iout,*) (L(j),j=1,3) return end subroutine inertia_tensor !----------------------------------------------------------------------------- subroutine angmom(cm,L) use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' real(kind=8) :: mscab real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp integer :: iti,inres,i,j,mnum ! Calculate the angular momentum do j=1,3 L(j)=0.0d0 enddo do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct-1 mnum=molnum(i) if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) enddo do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo call vecpr(pr(1),v(1),vp) do j=1,3 L(j)=L(j)+mp(mnum)*vp(j) ! print *,"HERE3",J,i,L(j),mp(mnum),Ip(mnum),mnum enddo do j=1,3 pr(j)=0.5d0*dc(j,i) pp(j)=0.5d0*d_t(j,i) enddo call vecpr(pr(1),pp(1),vp) ! print *,vp,"vp" do j=1,3 L(j)=L(j)+Ip(mnum)*vp(j) enddo enddo do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) inres=i+nres if (mnum.gt.4) then mscab=0.0d0 else mscab=msc(iti,mnum) endif do j=1,3 pr(j)=c(j,inres)-cm(j) enddo !endif if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo else do j=1,3 v(j)=incr(j) enddo endif ! print *,i,pr,"pr",v call vecpr(pr(1),v(1),vp) ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),& ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3) do j=1,3 L(j)=L(j)+mscab*vp(j) enddo ! write (iout,*) "L",(l(j),j=1,3) ! print *,"L",(l(j),j=1,3),i,vp(1) if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo call vecpr(dc(1,inres),d_t(1,inres),vp) do j=1,3 L(j)=L(j)+Isc(iti,mnum)*vp(j) enddo endif do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo enddo return end subroutine angmom !----------------------------------------------------------------------------- subroutine vcm_vel(vcm) use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.VAR' ! include 'COMMON.MD' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' real(kind=8),dimension(3) :: vcm,vv real(kind=8) :: summas,amas integer :: i,j,mnum do j=1,3 vcm(j)=0.0d0 vv(j)=d_t(j,0) enddo summas=0.0d0 do i=nnt,nct mnum=molnum(i) ! if (mnum.ge.4) mp(mnum)=msc(itype(i,mnum),mnum) if (i.lt.nct) then summas=summas+mp(mnum) do j=1,3 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i)) ! print *,i,j,vv(j),d_t(j,i) enddo endif ! if (mnum.ne.4) then amas=msc(iabs(itype(i,mnum)),mnum) ! else ! amas=0.0d0 ! endif summas=summas+amas if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then do j=1,3 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres)) enddo else do j=1,3 vcm(j)=vcm(j)+amas*vv(j) enddo endif do j=1,3 vv(j)=vv(j)+d_t(j,i) enddo enddo write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas do j=1,3 vcm(j)=vcm(j)/summas enddo return end subroutine vcm_vel #endif !----------------------------------------------------------------------------- ! rattle.F !----------------------------------------------------------------------------- subroutine rattle1 ! RATTLE algorithm for velocity Verlet - step 1, UNRES ! AL 9/24/04 use comm_przech use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef RATTLE ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' !el real(kind=8) :: gginv(2*nres,2*nres),& !el gdc(3,2*nres,2*nres) real(kind=8) :: dC_uncor(3,2*nres) !,& !el real(kind=8) :: Cmat(2*nres,2*nres) real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres !el common /przechowalnia/ GGinv,gdc,Cmat,nbond !el common /przechowalnia/ nbond integer :: max_rattle = 5 logical :: lprn = .false., lprn1 = .false., not_done real(kind=8) :: tol_rattle = 1.0d-5 integer :: ii,i,j,jj,l,ind,ind1,nres2 nres2=2*nres !el /common/ przechowalnia if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2)) if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2)) if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2)) !el-------- if (lprn) write (iout,*) "RATTLE1" nbond=nct-nnt do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) nbond=nbond+1 enddo ! Make a folded form of the Ginv-matrix ind=0 ii=0 do i=nnt,nct-1 ii=ii+1 do j=1,3 ind=ind+1 ind1=0 jj=0 do k=nnt,nct-1 jj=jj+1 do l=1,3 ind1=ind1+1 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1) enddo enddo do k=nnt,nct mnum=molnum(k) if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then jj=jj+1 do l=1,3 ind1=ind1+1 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1) enddo endif enddo enddo enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) ii=ii+1 do j=1,3 ind=ind+1 ind1=0 jj=0 do k=nnt,nct-1 jj=jj+1 do l=1,3 ind1=ind1+1 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1) enddo enddo do k=nnt,nct if (itype(k,1).ne.10) then jj=jj+1 do l=1,3 ind1=ind1+1 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1) enddo endif enddo enddo endif enddo if (lprn1) then write (iout,*) "Matrix GGinv" call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv) endif not_done=.true. iter=0 do while (not_done) iter=iter+1 if (iter.gt.max_rattle) then write (iout,*) "Error - too many iterations in RATTLE." stop endif ! Calculate the matrix C = GG**(-1) dC_old o dC ind1=0 do i=nnt,nct-1 ind1=ind1+1 do j=1,3 dC_uncor(j,ind1)=dC(j,i) enddo enddo do i=nnt,nct if (itype(i,1).ne.10) then ind1=ind1+1 do j=1,3 dC_uncor(j,ind1)=dC(j,i+nres) enddo endif enddo do i=1,nbond ind=0 do k=nnt,nct-1 ind=ind+1 do j=1,3 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k) enddo enddo do k=nnt,nct if (itype(k,1).ne.10) then ind=ind+1 do j=1,3 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres) enddo endif enddo enddo ! Calculate deviations from standard virtual-bond lengths ind=0 do i=nnt,nct-1 ind=ind+1 x(ind)=vbld(i+1)**2-vbl**2 enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2 endif enddo if (lprn) then write (iout,*) "Coordinates and violations" do i=1,nbond write(iout,'(i5,3f10.5,5x,e15.5)') & i,(dC_uncor(j,i),j=1,3),x(i) enddo write (iout,*) "Velocities and violations" ind=0 do i=nnt,nct-1 ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i)) enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t_new(j,i+nres),j=1,3),& scalar(d_t_new(1,i+nres),dC_old(1,i+nres)) endif enddo ! write (iout,*) "gdc" ! do i=1,nbond ! write (iout,*) "i",i ! do j=1,nbond ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3) ! enddo ! enddo endif xmax=dabs(x(1)) do i=2,nbond if (dabs(x(i)).gt.xmax) then xmax=dabs(x(i)) endif enddo if (xmax.lt.tol_rattle) then not_done=.false. goto 100 endif ! Calculate the matrix of the system of equations do i=1,nbond do j=1,nbond Cmat(i,j)=0.0d0 do k=1,3 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j) enddo enddo enddo if (lprn1) then write (iout,*) "Matrix Cmat" call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat) endif call gauss(Cmat,X,MAXRES2,nbond,1,*10) ! Add constraint term to positions ind=0 do i=nnt,nct-1 ind=ind+1 do j=1,3 xx=0.0d0 do ii=1,nbond xx = xx+x(ii)*gdc(j,ind,ii) enddo xx=0.5d0*xx dC(j,i)=dC(j,i)-xx d_t_new(j,i)=d_t_new(j,i)-xx/d_time enddo enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 do j=1,3 xx=0.0d0 do ii=1,nbond xx = xx+x(ii)*gdc(j,ind,ii) enddo xx=0.5d0*xx dC(j,i+nres)=dC(j,i+nres)-xx d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time enddo endif enddo ! Rebuild the chain using the new coordinates call chainbuild_cart if (lprn) then write (iout,*) "New coordinates, Lagrange multipliers,",& " and differences between actual and standard bond lengths" ind=0 do i=nnt,nct-1 ind=ind+1 xx=vbld(i+1)**2-vbl**2 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i),j=1,3),x(ind),xx enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) ind=ind+1 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i+nres),j=1,3),x(ind),xx endif enddo write (iout,*) "Velocities and violations" ind=0 do i=nnt,nct-1 ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i)) enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t_new(j,i+nres),j=1,3),& scalar(d_t_new(1,i+nres),dC_old(1,i+nres)) endif enddo endif enddo 100 continue return 10 write (iout,*) "Error - singularity in solving the system",& " of equations for Lagrange multipliers." stop #else write (iout,*) & "RATTLE inactive; use -DRATTLE switch at compile time." stop #endif end subroutine rattle1 !----------------------------------------------------------------------------- subroutine rattle2 ! RATTLE algorithm for velocity Verlet - step 2, UNRES ! AL 9/24/04 use comm_przech use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef RATTLE ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' !el real(kind=8) :: gginv(2*nres,2*nres),& !el gdc(3,2*nres,2*nres) real(kind=8) :: dC_uncor(3,2*nres) !,& !el Cmat(2*nres,2*nres) real(kind=8) :: x(2*nres) !maxres2=2*maxres !el common /przechowalnia/ GGinv,gdc,Cmat,nbond !el common /przechowalnia/ nbond integer :: max_rattle = 5 logical :: lprn = .false., lprn1 = .false., not_done real(kind=8) :: tol_rattle = 1.0d-5 integer :: nres2 nres2=2*nres !el /common/ przechowalnia if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2)) if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2)) if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2)) !el-------- if (lprn) write (iout,*) "RATTLE2" if (lprn) write (iout,*) "Velocity correction" ! Calculate the matrix G dC do i=1,nbond ind=0 do k=nnt,nct-1 ind=ind+1 do j=1,3 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k) enddo enddo do k=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then ind=ind+1 do j=1,3 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres) enddo endif enddo enddo ! if (lprn) then ! write (iout,*) "gdc" ! do i=1,nbond ! write (iout,*) "i",i ! do j=1,nbond ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3) ! enddo ! enddo ! endif ! Calculate the matrix of the system of equations ind=0 do i=nnt,nct-1 ind=ind+1 do j=1,nbond Cmat(ind,j)=0.0d0 do k=1,3 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j) enddo enddo enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then ind=ind+1 do j=1,nbond Cmat(ind,j)=0.0d0 do k=1,3 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j) enddo enddo endif enddo ! Calculate the scalar product dC o d_t_new ind=0 do i=nnt,nct-1 ind=ind+1 x(ind)=scalar(d_t(1,i),dC(1,i)) enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then ind=ind+1 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres)) endif enddo if (lprn) then write (iout,*) "Velocities and violations" ind=0 do i=nnt,nct-1 ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i,ind,(d_t(j,i),j=1,3),x(ind) enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind) endif enddo endif xmax=dabs(x(1)) do i=2,nbond if (dabs(x(i)).gt.xmax) then xmax=dabs(x(i)) endif enddo if (xmax.lt.tol_rattle) then not_done=.false. goto 100 endif if (lprn1) then write (iout,*) "Matrix Cmat" call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat) endif call gauss(Cmat,X,MAXRES2,nbond,1,*10) ! Add constraint term to velocities ind=0 do i=nnt,nct-1 ind=ind+1 do j=1,3 xx=0.0d0 do ii=1,nbond xx = xx+x(ii)*gdc(j,ind,ii) enddo d_t(j,i)=d_t(j,i)-xx enddo enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) then ind=ind+1 do j=1,3 xx=0.0d0 do ii=1,nbond xx = xx+x(ii)*gdc(j,ind,ii) enddo d_t(j,i+nres)=d_t(j,i+nres)-xx enddo endif enddo if (lprn) then write (iout,*) & "New velocities, Lagrange multipliers violations" ind=0 do i=nnt,nct-1 ind=ind+1 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') & i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i)) enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& .and.(mnum.lt.4)) ind=ind+1 write (iout,'(2i5,3f10.5,5x,2e15.5)') & i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),& scalar(d_t(1,i+nres),dC(1,i+nres)) endif enddo endif 100 continue return 10 write (iout,*) "Error - singularity in solving the system",& " of equations for Lagrange multipliers." stop #else write (iout,*) & "RATTLE inactive; use -DRATTLE option at compile time." stop #endif end subroutine rattle2 !----------------------------------------------------------------------------- subroutine rattle_brown ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES ! AL 9/24/04 use comm_przech use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef RATTLE ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' !el real(kind=8) :: gginv(2*nres,2*nres),& !el gdc(3,2*nres,2*nres) real(kind=8) :: dC_uncor(3,2*nres) !,& !el real(kind=8) :: Cmat(2*nres,2*nres) real(kind=8) :: x(2*nres) !maxres2=2*maxres !el common /przechowalnia/ GGinv,gdc,Cmat,nbond !el common /przechowalnia/ nbond integer :: max_rattle = 5 logical :: lprn = .false., lprn1 = .false., not_done real(kind=8) :: tol_rattle = 1.0d-5 integer :: nres2 nres2=2*nres !el /common/ przechowalnia if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2)) if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2)) if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2)) !el-------- if (lprn) write (iout,*) "RATTLE_BROWN" nbond=nct-nnt do i=nnt,nct if (itype(i,1).ne.10) nbond=nbond+1 enddo ! Make a folded form of the Ginv-matrix ind=0 ii=0 do i=nnt,nct-1 ii=ii+1 do j=1,3 ind=ind+1 ind1=0 jj=0 do k=nnt,nct-1 jj=jj+1 do l=1,3 ind1=ind1+1 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1) enddo enddo do k=nnt,nct if (itype(k,1).ne.10) then jj=jj+1 do l=1,3 ind1=ind1+1 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1) enddo endif enddo enddo enddo do i=nnt,nct if (itype(i,1).ne.10) then ii=ii+1 do j=1,3 ind=ind+1 ind1=0 jj=0 do k=nnt,nct-1 jj=jj+1 do l=1,3 ind1=ind1+1 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1) enddo enddo do k=nnt,nct if (itype(k,1).ne.10) then jj=jj+1 do l=1,3 ind1=ind1+1 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1) enddo endif enddo enddo endif enddo if (lprn1) then write (iout,*) "Matrix GGinv" call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv) endif not_done=.true. iter=0 do while (not_done) iter=iter+1 if (iter.gt.max_rattle) then write (iout,*) "Error - too many iterations in RATTLE." stop endif ! Calculate the matrix C = GG**(-1) dC_old o dC ind1=0 do i=nnt,nct-1 ind1=ind1+1 do j=1,3 dC_uncor(j,ind1)=dC(j,i) enddo enddo do i=nnt,nct if (itype(i,1).ne.10) then ind1=ind1+1 do j=1,3 dC_uncor(j,ind1)=dC(j,i+nres) enddo endif enddo do i=1,nbond ind=0 do k=nnt,nct-1 ind=ind+1 do j=1,3 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k) enddo enddo do k=nnt,nct if (itype(k,1).ne.10) then ind=ind+1 do j=1,3 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres) enddo endif enddo enddo ! Calculate deviations from standard virtual-bond lengths ind=0 do i=nnt,nct-1 ind=ind+1 x(ind)=vbld(i+1)**2-vbl**2 enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2 endif enddo if (lprn) then write (iout,*) "Coordinates and violations" do i=1,nbond write(iout,'(i5,3f10.5,5x,e15.5)') & i,(dC_uncor(j,i),j=1,3),x(i) enddo write (iout,*) "Velocities and violations" ind=0 do i=nnt,nct-1 ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i)) enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t(j,i+nres),j=1,3),& scalar(d_t(1,i+nres),dC_old(1,i+nres)) endif enddo write (iout,*) "gdc" do i=1,nbond write (iout,*) "i",i do j=1,nbond write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3) enddo enddo endif xmax=dabs(x(1)) do i=2,nbond if (dabs(x(i)).gt.xmax) then xmax=dabs(x(i)) endif enddo if (xmax.lt.tol_rattle) then not_done=.false. goto 100 endif ! Calculate the matrix of the system of equations do i=1,nbond do j=1,nbond Cmat(i,j)=0.0d0 do k=1,3 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j) enddo enddo enddo if (lprn1) then write (iout,*) "Matrix Cmat" call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat) endif call gauss(Cmat,X,MAXRES2,nbond,1,*10) ! Add constraint term to positions ind=0 do i=nnt,nct-1 ind=ind+1 do j=1,3 xx=0.0d0 do ii=1,nbond xx = xx+x(ii)*gdc(j,ind,ii) enddo xx=-0.5d0*xx d_t(j,i)=d_t(j,i)+xx/d_time dC(j,i)=dC(j,i)+xx enddo enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 do j=1,3 xx=0.0d0 do ii=1,nbond xx = xx+x(ii)*gdc(j,ind,ii) enddo xx=-0.5d0*xx d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time dC(j,i+nres)=dC(j,i+nres)+xx enddo endif enddo ! Rebuild the chain using the new coordinates call chainbuild_cart if (lprn) then write (iout,*) "New coordinates, Lagrange multipliers,",& " and differences between actual and standard bond lengths" ind=0 do i=nnt,nct-1 ind=ind+1 xx=vbld(i+1)**2-vbl**2 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i),j=1,3),x(ind),xx enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i+nres),j=1,3),x(ind),xx endif enddo write (iout,*) "Velocities and violations" ind=0 do i=nnt,nct-1 ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i)) enddo do i=nnt,nct if (itype(i,1).ne.10) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t_new(j,i+nres),j=1,3),& scalar(d_t_new(1,i+nres),dC_old(1,i+nres)) endif enddo endif enddo 100 continue return 10 write (iout,*) "Error - singularity in solving the system",& " of equations for Lagrange multipliers." stop #else write (iout,*) & "RATTLE inactive; use -DRATTLE option at compile time" stop #endif end subroutine rattle_brown !----------------------------------------------------------------------------- ! stochfric.F !----------------------------------------------------------------------------- subroutine friction_force use energy_data use REMD_data use comm_syfek ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.IOUNITS' !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres !el common /syfek/ gamvec #ifdef FIVEDIAG integer iposc,ichain,n,innt,inct double precision rs(nres*2) real(kind=8) ::v_work(3,6*nres),vvec(2*nres) #else real(kind=8) :: v_work(6*nres) #endif real(kind=8) :: vv(3),vvtot(3,nres)!,v_work(6*nres) !,& !el ginvfric(2*nres,2*nres) !maxres2=2*maxres !el common /przechowalnia/ ginvfric logical :: lprn, checkmode integer :: i,j,ind,k,nres2,nres6,mnum nres2=2*nres nres6=6*nres lprn=.false. checkmode=.false. ! if (large) lprn=.true. ! if (large) checkmode=.true. #ifdef FIVEDIAG !c Here accelerations due to friction forces are computed right after forces. d_t_work(:6*nres)=0.0d0 do j=1,3 v_work(j,1)=d_t(j,0) v_work(j,nnt)=d_t(j,0) enddo do i=nnt+1,nct do j=1,3 v_work(j,i)=v_work(j,i-1)+d_t(j,i-1) enddo enddo do i=nnt,nct mnum=molnum(i) if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then do j=1,3 v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres) enddo endif enddo #ifdef DEBUG write (iout,*) "v_work" do i=1,2*nres write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3) enddo #endif do j=1,3 ind=0 do ichain=1,nchain n=dimen_chain(ichain) iposc=iposd_chain(ichain) !c write (iout,*) "friction_force j",j," ichain",ichain, !c & " n",n," iposc",iposc,iposc+n-1 innt=chain_border(1,ichain) inct=chain_border(2,ichain) !c diagnostics !c innt=chain_border(1,1) !c inct=chain_border(2,1) do i=innt,inct mnum=molnum(i) vvec(ind+1)=v_work(j,i) ind=ind+1 ! if (iabs(itype(i)).ne.10) then if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then vvec(ind+1)=v_work(j,i+nres) ind=ind+1 endif enddo #ifdef DEBUG write (iout,*) "vvec ind",ind," n",n write (iout,'(f10.5)') (vvec(i),i=iposc,ind) #endif !c write (iout,*) "chain",i," ind",ind," n",n #ifdef TIMING #ifdef MPI time01=MPI_Wtime() #else time01=tcpu() #endif #endif ! if (large) print *,"before fivediagmult" call fivediagmult(n,DMfric(iposc),DU1fric(iposc),& DU2fric(iposc),vvec(iposc),rs) ! if (large) print *,"after fivediagmult" #ifdef TIMING #ifdef MPI time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01 #else time_fricmatmult=time_fricmatmult+tcpu()-time01 #endif #endif #ifdef DEBUG write (iout,*) "rs" write (iout,'(f10.5)') (rs(i),i=1,n) #endif do i=iposc,iposc+n-1 ! write (iout,*) "ichain",ichain," i",i," j",j,& ! "index",3*(i-1)+j,"rs",rs(i-iposc+1) fric_work(3*(i-1)+j)=-rs(i-iposc+1) enddo enddo enddo #ifdef DEBUG write (iout,*) "Vector fric_work dimen3",dimen3 write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3) #endif #else if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6) if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres do i=0,nres2 do j=1,3 friction(j,i)=0.0d0 enddo enddo do j=1,3 d_t_work(j)=d_t(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 d_t_work(ind+j)=d_t(j,i) enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& .and.(mnum.lt.4)) then do j=1,3 d_t_work(ind+j)=d_t(j,i+nres) enddo ind=ind+3 endif enddo call fricmat_mult(d_t_work,fric_work) if (.not.checkmode) return if (lprn) then write (iout,*) "d_t_work and fric_work" do i=1,3*dimen write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i) enddo endif do j=1,3 friction(j,0)=fric_work(j) enddo ind=3 do i=nnt,nct-1 do j=1,3 friction(j,i)=fric_work(ind+j) enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& .and.(mnum.lt.4)) then ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then do j=1,3 friction(j,i+nres)=fric_work(ind+j) enddo ind=ind+3 endif enddo if (lprn) then write(iout,*) "Friction backbone" do i=0,nct-1 write(iout,'(i5,3e15.5,5x,3e15.5)') & i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3) enddo write(iout,*) "Friction side chain" do i=nnt,nct write(iout,'(i5,3e15.5,5x,3e15.5)') & i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3) enddo endif if (lprn) then do j=1,3 vv(j)=d_t(j,0) enddo do i=nnt,nct do j=1,3 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i) vvtot(j,i+nres)=vv(j)+d_t(j,i+nres) vv(j)=vv(j)+d_t(j,i) enddo enddo write (iout,*) "vvtot backbone and sidechain" do i=nnt,nct write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),& (vvtot(j,i+nres),j=1,3) enddo ind=0 do i=nnt,nct-1 do j=1,3 v_work(ind+j)=vvtot(j,i) enddo ind=ind+3 enddo do i=nnt,nct do j=1,3 v_work(ind+j)=vvtot(j,i+nres) enddo ind=ind+3 enddo write (iout,*) "v_work gamvec and site-based friction forces" do i=1,dimen1 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),& gamvec(i)*v_work(i) enddo ! do i=1,dimen ! fric_work1(i)=0.0d0 ! do j=1,dimen1 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j) ! enddo ! enddo ! write (iout,*) "fric_work and fric_work1" ! do i=1,dimen ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i) ! enddo do i=1,dimen do j=1,dimen ginvfric(i,j)=0.0d0 do k=1,dimen ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j) enddo enddo enddo write (iout,*) "ginvfric" do i=1,dimen write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen) enddo write (iout,*) "symmetry check" do i=1,dimen do j=1,i-1 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i) enddo enddo endif #endif return end subroutine friction_force !----------------------------------------------------------------------------- subroutine setup_fricmat ! use MPI use energy_data use control_data, only:time_Bcast use control, only:tcpu use comm_syfek ! implicit real*8 (a-h,o-z) #ifdef MPI use MPI_data include 'mpif.h' real(kind=8) :: time00 #endif ! include 'DIMENSIONS' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.MD' ! include 'COMMON.SETUP' ! include 'COMMON.TIME1' ! integer licznik /0/ ! save licznik !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.IOUNITS' integer :: IERROR integer :: i,j,ind,ind1,m,ichain,innt,inct logical :: lprn = .true. real(kind=8) :: dtdi !el ,gamvec(2*nres) !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy ! real(kind=8),allocatable,dimension(:,:) :: fcopy !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2)) !el common /syfek/ gamvec real(kind=8) :: work(8*2*nres) integer :: iwork(2*nres) !el common /przechowalnia/ ginvfric,Ghalf,fcopy integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum nres2=2*nres nres6=6*nres #ifdef MPI #ifndef FIVEDIAG if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres if (fg_rank.ne.king) goto 10 #endif #endif ! nres2=2*nres ! nres6=6*nres if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2) #ifndef FIVEDIAG if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) #endif #ifdef FIVEDIAG if (.not.allocated(DMfric)) allocate(DMfric(nres2)) if (.not.allocated(DU1fric)) allocate(DU1fric(nres2)) if (.not.allocated(DU2fric)) allocate(DU2fric(nres2)) ! Load the friction coefficients corresponding to peptide groups ind1=0 do i=nnt,nct-1 mnum=molnum(i) ind1=ind1+1 gamvec(ind1)=gamp(mnum) enddo !HERE TEST if (molnum(nct).eq.5) then mnum=molnum(i) ind1=ind1+1 gamvec(ind1)=gamp(mnum) endif ! Load the friction coefficients corresponding to side chains m=nct-nnt ind=0 do j=1,2 gamsc(ntyp1_molec(j),j)=1.0d0 enddo do i=nnt,nct mnum=molnum(i) ind=ind+1 ii = ind+m iti=itype(i,mnum) gamvec(ii)=gamsc(iabs(iti),mnum) enddo if (surfarea) call sdarea(gamvec) DMfric=0.0d0 DU1fric=0.0d0 DU2fric=0.0d0 ind=1 do ichain=1,nchain innt=chain_border(1,ichain) inct=chain_border(2,ichain) !c write (iout,*) "ichain",ichain," innt",innt," inct",inct !c DMfric part mnum=molnum(innt) DMfric(ind)=gamvec(innt-nnt+1)/4 if (iabs(itype(innt,1)).eq.10.or.mnum.gt.2) then DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1) ind=ind+1 else DMfric(ind+1)=gamvec(m+innt-nnt+1) ind=ind+2 endif !c write (iout,*) "DMfric init ind",ind !c DMfric do i=innt+1,inct-1 mnum=molnum(i) DMfric(ind)=gamvec(i-nnt+1)/2 if (iabs(itype(i,1)).eq.10.or.mnum.gt.2) then DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1) ind=ind+1 else DMfric(ind+1)=gamvec(m+i-nnt+1) ind=ind+2 endif enddo !c write (iout,*) "DMfric endloop ind",ind if (inct.gt.innt) then DMfric(ind)=gamvec(inct-1-nnt+1)/4 mnum=molnum(inct) if (iabs(itype(inct,1)).eq.10.or.mnum.gt.2) then DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1) ind=ind+1 else DMfric(ind+1)=gamvec(inct+m-nnt+1) ind=ind+2 endif endif !c write (iout,*) "DMfric end ind",ind enddo !c DU1fric part do ichain=1,nchain ind=iposd_chain(ichain) innt=chain_border(1,ichain) inct=chain_border(2,ichain) do i=innt,inct mnum=molnum(i) if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then ind=ind+2 else DU1fric(ind)=gamvec(i-nnt+1)/4 ind=ind+1 endif enddo enddo !c DU2fric part do ichain=1,nchain ind=iposd_chain(ichain) innt=chain_border(1,ichain) inct=chain_border(2,ichain) do i=innt,inct-1 mnum=molnum(i) if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then DU2fric(ind)=gamvec(i-nnt+1)/4 DU2fric(ind+1)=0.0d0 ind=ind+2 else DU2fric(ind)=0.0d0 ind=ind+1 endif enddo enddo if (lprn) then write(iout,*)"The upper part of the five-diagonal friction matrix" do ichain=1,nchain write (iout,'(a,i5)') 'Chain',ichain innt=iposd_chain(ichain) inct=iposd_chain(ichain)+dimen_chain(ichain)-1 do i=innt,inct if (i.lt.inct-1) then write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),& DU2fric(i) else if (i.eq.inct-1) then write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i) else write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i) endif enddo enddo endif 10 continue #else ! Zeroing out fricmat do i=1,dimen do j=1,dimen fricmat(i,j)=0.0d0 enddo enddo ! Load the friction coefficients corresponding to peptide groups ind1=0 do i=nnt,nct-1 mnum=molnum(i) ind1=ind1+1 gamvec(ind1)=gamp(mnum) enddo !HERE TEST if (molnum(nct).eq.5) then mnum=molnum(i) ind1=ind1+1 gamvec(ind1)=gamp(mnum) endif ! Load the friction coefficients corresponding to side chains m=nct-nnt ind=0 do j=1,2 gamsc(ntyp1_molec(j),j)=1.0d0 enddo do i=nnt,nct mnum=molnum(i) ind=ind+1 ii = ind+m iti=itype(i,mnum) gamvec(ii)=gamsc(iabs(iti),mnum) enddo if (surfarea) call sdarea(gamvec) ! if (lprn) then ! write (iout,*) "Matrix A and vector gamma" ! do i=1,dimen1 ! write (iout,'(i2,$)') i ! do j=1,dimen ! write (iout,'(f4.1,$)') A(i,j) ! enddo ! write (iout,'(f8.3)') gamvec(i) ! enddo ! endif if (lprn) then write (iout,*) "Vector gamvec" do i=1,dimen1 write (iout,'(i5,f10.5)') i, gamvec(i) enddo endif ! The friction matrix do k=1,dimen do i=1,dimen dtdi=0.0d0 do j=1,dimen1 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j) enddo fricmat(k,i)=dtdi enddo enddo if (lprn) then write (iout,'(//a)') "Matrix fricmat" call matout2(dimen,dimen,nres2,nres2,fricmat) endif if (lang.eq.2 .or. lang.eq.3) then ! Mass-scale the friction matrix if non-direct integration will be performed do i=1,dimen do j=1,dimen Ginvfric(i,j)=0.0d0 do k=1,dimen do l=1,dimen Ginvfric(i,j)=Ginvfric(i,j)+ & Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l) enddo enddo enddo enddo ! Diagonalize the friction matrix ind=0 do i=1,dimen do j=1,i ind=ind+1 Ghalf(ind)=Ginvfric(i,j) enddo enddo call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,& ierr,iwork) if (lprn) then write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",& " mass-scaled friction matrix" call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam) endif ! Precompute matrices for tinker stochastic integrator #ifndef LANG0 do i=1,dimen do j=1,dimen mt1(i,j)=0.0d0 mt2(i,j)=0.0d0 do k=1,dimen mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j) mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j) enddo mt3(j,i)=mt1(i,j) enddo enddo #endif else if (lang.eq.4) then ! Diagonalize the friction matrix ind=0 do i=1,dimen do j=1,i ind=ind+1 Ghalf(ind)=fricmat(i,j) enddo enddo call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,& ierr,iwork) if (lprn) then write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",& " friction matrix" call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam) endif ! Determine the number of zero eigenvalues of the friction matrix nzero=max0(dimen-dimen1,0) ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen) ! nzero=nzero+1 ! enddo write (iout,*) "Number of zero eigenvalues:",nzero do i=1,dimen do j=1,dimen fricmat(i,j)=0.0d0 do k=nzero+1,dimen fricmat(i,j)=fricmat(i,j) & +fricvec(i,k)*fricvec(j,k)/fricgam(k) enddo enddo enddo if (lprn) then write (iout,'(//a)') "Generalized inverse of fricmat" call matout(dimen,dimen,nres6,nres6,fricmat) endif endif #ifdef MPI 10 continue if (nfgtasks.gt.1) then if (fg_rank.eq.0) then ! The matching BROADCAST for fg processors is called in ERGASTULUM #ifdef MPI time00=MPI_Wtime() #else time00=tcpu() #endif call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR) #ifdef MPI time_Bcast=time_Bcast+MPI_Wtime()-time00 #else time_Bcast=time_Bcast+tcpu()-time00 #endif ! print *,"Processor",myrank, ! & " BROADCAST iorder in SETUP_FRICMAT" endif ! licznik=licznik+1 write (iout,*) "setup_fricmat licznik"!,licznik !sp #ifdef MPI time00=MPI_Wtime() #else time00=tcpu() #endif ! Scatter the friction matrix call MPI_Scatterv(fricmat(1,1),nginv_counts(0),& nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),& myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR) #ifdef TIMING #ifdef MPI time_scatter=time_scatter+MPI_Wtime()-time00 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00 #else time_scatter=time_scatter+tcpu()-time00 time_scatter_fmat=time_scatter_fmat+tcpu()-time00 #endif #endif do i=1,dimen do j=1,2*my_ng_count fricmat(j,i)=fcopy(i,j) enddo enddo ! write (iout,*) "My chunk of fricmat" ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy) endif #endif #endif return end subroutine setup_fricmat !----------------------------------------------------------------------------- subroutine stochastic_force(stochforcvec) use energy_data use random, only:anorm_distr ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' use control, only: tcpu use control_data #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.MD' ! include 'COMMON.TIME1' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.IOUNITS' real(kind=8) :: x,sig,lowb,highb real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres) real(kind=8) :: time00 logical :: lprn = .false. integer :: i,j,ind,mnum #ifdef FIVEDIAG integer ichain,innt,inct,iposc #endif do i=0,2*nres do j=1,3 stochforc(j,i)=0.0d0 enddo enddo x=0.0d0 #ifdef MPI time00=MPI_Wtime() #else time00=tcpu() #endif ! Compute the stochastic forces acting on bodies. Store in force. do i=nnt,nct-1 sig=stdforcp(i) lowb=-5*sig highb=5*sig do j=1,3 force(j,i)=anorm_distr(x,sig,lowb,highb) enddo enddo do i=nnt,nct sig2=stdforcsc(i) lowb2=-5*sig2 highb2=5*sig2 do j=1,3 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2) enddo enddo #ifdef MPI time_fsample=time_fsample+MPI_Wtime()-time00 #else time_fsample=time_fsample+tcpu()-time00 #endif #ifdef FIVEDIAG ind=0 do ichain=1,nchain innt=chain_border(1,ichain) inct=chain_border(2,ichain) iposc=iposd_chain(ichain) !c for debugging only !c innt=chain_border(1,1) !c inct=chain_border(2,1) !c iposc=iposd_chain(1) !c write (iout,*)"stochastic_force ichain=",ichain," innt",innt, !c & " inct",inct," iposc",iposc do j=1,3 stochforcvec(ind+j)=0.5d0*force(j,innt) enddo if (iabs(itype(innt,molnum(innt))).eq.10.or.molnum(innt).ge.3) then do j=1,3 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres) enddo ind=ind+3 else ind=ind+3 do j=1,3 stochforcvec(ind+j)=force(j,innt+nres) enddo ind=ind+3 endif do i=innt+1,inct-1 do j=1,3 stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1)) enddo if (iabs(itype(i,1).eq.10).or.molnum(i).ge.3) then do j=1,3 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres) enddo ind=ind+3 else ind=ind+3 do j=1,3 stochforcvec(ind+j)=force(j,i+nres) enddo ind=ind+3 endif enddo do j=1,3 stochforcvec(ind+j)=0.5d0*force(j,inct-1) enddo if (iabs(itype(inct,molnum(inct))).eq.10.or.molnum(inct).ge.3) then do j=1,3 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres) enddo ind=ind+3 else ind=ind+3 do j=1,3 stochforcvec(ind+j)=force(j,inct+nres) enddo ind=ind+3 endif !c write (iout,*) "chain",ichain," ind",ind enddo #ifdef DEBUG write (iout,*) "stochforcvec" write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind) #endif #else ! Compute the stochastic forces acting on virtual-bond vectors. do j=1,3 ff(j)=0.0d0 enddo do i=nct-1,nnt,-1 do j=1,3 stochforc(j,i)=ff(j)+0.5d0*force(j,i) enddo do j=1,3 ff(j)=ff(j)+force(j,i) enddo ! if (itype(i+1,1).ne.ntyp1) then mnum=molnum(i) if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then do j=1,3 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1) ff(j)=ff(j)+force(j,i+nres+1) enddo endif enddo do j=1,3 stochforc(j,0)=ff(j)+force(j,nnt+nres) enddo do i=nnt,nct mnum=molnum(i) if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& .and.(mnum.lt.4)) then ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then do j=1,3 stochforc(j,i+nres)=force(j,i+nres) enddo endif enddo do j=1,3 stochforcvec(j)=stochforc(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 stochforcvec(ind+j)=stochforc(j,i) enddo ind=ind+3 enddo do i=nnt,nct mnum=molnum(i) if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& .and.(mnum.lt.4)) then ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then do j=1,3 stochforcvec(ind+j)=stochforc(j,i+nres) enddo ind=ind+3 endif enddo if (lprn) then write (iout,*) "stochforcvec" do i=1,3*dimen write(iout,'(i5,e15.5)') i,stochforcvec(i) enddo write(iout,*) "Stochastic forces backbone" do i=0,nct-1 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3) enddo write(iout,*) "Stochastic forces side chain" do i=nnt,nct write(iout,'(i5,3e15.5)') & i,(stochforc(j,i+nres),j=1,3) enddo endif if (lprn) then ind=0 do i=nnt,nct-1 write (iout,*) i,ind do j=1,3 forcvec(ind+j)=force(j,i) enddo ind=ind+3 enddo do i=nnt,nct write (iout,*) i,ind do j=1,3 forcvec(j+ind)=force(j,i+nres) enddo ind=ind+3 enddo write (iout,*) "forcvec" ind=0 do i=nnt,nct-1 do j=1,3 write (iout,'(2i3,2f10.5)') i,j,force(j,i),& forcvec(ind+j) enddo ind=ind+3 enddo do i=nnt,nct do j=1,3 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),& forcvec(ind+j) enddo ind=ind+3 enddo endif #endif return end subroutine stochastic_force !----------------------------------------------------------------------------- subroutine sdarea(gamvec) ! ! Scale the friction coefficients according to solvent accessible surface areas ! Code adapted from TINKER ! AL 9/3/04 ! use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' ! include 'COMMON.MD' !#ifndef LANG0 ! include 'COMMON.LANGEVIN' !#else ! include 'COMMON.LANGEVIN.lang0' !#endif ! include 'COMMON.CHAIN' ! include 'COMMON.DERIV' ! include 'COMMON.GEO' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2) real(kind=8),parameter :: twosix = 1.122462048309372981d0 logical :: lprn = .false. real(kind=8) :: probe,area,ratio integer :: i,j,ind,iti,mnum ! ! determine new friction coefficients every few SD steps ! ! set the atomic radii to estimates of sigma values ! ! print *,"Entered sdarea" probe = 0.0d0 do i=1,2*nres radius(i)=0.0d0 enddo ! Load peptide group radii do i=nnt,nct-1 mnum=molnum(i) radius(i)=pstok(mnum) enddo ! Load side chain radii do i=nnt,nct mnum=molnum(i) iti=itype(i,mnum) radius(i+nres)=restok(iti,mnum) enddo ! do i=1,2*nres ! write (iout,*) "i",i," radius",radius(i) ! enddo do i = 1, 2*nres radius(i) = radius(i) / twosix if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe end do ! ! scale atomic friction coefficients by accessible area ! if (lprn) write (iout,*) & "Original gammas, surface areas, scaling factors, new gammas, ",& "std's of stochastic forces" ind=0 do i=nnt,nct-1 if (radius(i).gt.0.0d0) then call surfatom (i,area,radius) ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1) if (lprn) write (iout,'(i5,3f10.5,$)') & i,gamvec(ind+1),area,ratio do j=1,3 ind=ind+1 gamvec(ind) = ratio * gamvec(ind) enddo mnum=molnum(i) stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind)) if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i) endif enddo do i=nnt,nct if (radius(i+nres).gt.0.0d0) then call surfatom (i+nres,area,radius) ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1) if (lprn) write (iout,'(i5,3f10.5,$)') & i,gamvec(ind+1),area,ratio do j=1,3 ind=ind+1 gamvec(ind) = ratio * gamvec(ind) enddo mnum=molnum(i) stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind)) if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i) endif enddo return end subroutine sdarea !----------------------------------------------------------------------------- ! surfatom.f !----------------------------------------------------------------------------- ! ! ! ################################################### ! ## COPYRIGHT (C) 1996 by Jay William Ponder ## ! ## All Rights Reserved ## ! ################################################### ! ! ################################################################ ! ## ## ! ## subroutine surfatom -- exposed surface area of an atom ## ! ## ## ! ################################################################ ! ! ! "surfatom" performs an analytical computation of the surface ! area of a specified atom; a simplified version of "surface" ! ! literature references: ! ! T. J. Richmond, "Solvent Accessible Surface Area and ! Excluded Volume in Proteins", Journal of Molecular Biology, ! 178, 63-89 (1984) ! ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters ! Applied to Molecular Dynamics of Proteins in Solution", ! Protein Science, 1, 227-235 (1992) ! ! variables and parameters: ! ! ir number of atom for which area is desired ! area accessible surface area of the atom ! radius radii of each of the individual atoms ! ! subroutine surfatom(ir,area,radius) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'sizes.i' ! include 'COMMON.GEO' ! include 'COMMON.IOUNITS' ! integer :: nres, integer :: nsup,nstart_sup ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2), ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2), ! & dc_work(MAXRES6),nres,nres0 integer,parameter :: maxarc=300 integer :: i,j,k,m integer :: ii,ib,jb integer :: io,ir integer :: mi,ni,narc integer :: key(maxarc) integer :: intag(maxarc) integer :: intag1(maxarc) real(kind=8) :: area,arcsum real(kind=8) :: arclen,exang real(kind=8) :: delta,delta2 real(kind=8) :: eps,rmove real(kind=8) :: xr,yr,zr real(kind=8) :: rr,rrsq real(kind=8) :: rplus,rminus real(kind=8) :: axx,axy,axz real(kind=8) :: ayx,ayy real(kind=8) :: azx,azy,azz real(kind=8) :: uxj,uyj,uzj real(kind=8) :: tx,ty,tz real(kind=8) :: txb,tyb,td real(kind=8) :: tr2,tr,txr,tyr real(kind=8) :: tk1,tk2 real(kind=8) :: thec,the,t,tb real(kind=8) :: txk,tyk,tzk real(kind=8) :: t1,ti,tf,tt real(kind=8) :: txj,tyj,tzj real(kind=8) :: ccsq,cc,xysq real(kind=8) :: bsqk,bk,cosine real(kind=8) :: dsqj,gi,pix2 real(kind=8) :: therk,dk,gk real(kind=8) :: risqk,rik real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2) real(kind=8) :: ri(maxarc),risq(maxarc) real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc) real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc) real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc) real(kind=8) :: dsq(maxarc),bsq(maxarc) real(kind=8) :: dsq1(maxarc),bsq1(maxarc) real(kind=8) :: arci(maxarc),arcf(maxarc) real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc) real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc) real(kind=8) :: kent(maxarc),kout(maxarc) real(kind=8) :: ther(maxarc) logical :: moved,top logical :: omit(maxarc) ! ! include 'sizes.i' ! maxatm = 2*nres !maxres2 maxres2=2*maxres ! maxlight = 8*maxatm ! maxbnd = 2*maxatm ! maxang = 3*maxatm ! maxtors = 4*maxatm ! ! zero out the surface area for the sphere of interest ! area = 0.0d0 ! write (2,*) "ir",ir," radius",radius(ir) if (radius(ir) .eq. 0.0d0) return ! ! set the overlap significance and connectivity shift ! pix2 = 2.0d0 * pi delta = 1.0d-8 delta2 = delta * delta eps = 1.0d-8 moved = .false. rmove = 1.0d-8 ! ! store coordinates and radius of the sphere of interest ! xr = c(1,ir) yr = c(2,ir) zr = c(3,ir) rr = radius(ir) rrsq = rr * rr ! ! initialize values of some counters and summations ! 10 continue io = 0 jb = 0 ib = 0 arclen = 0.0d0 exang = 0.0d0 ! ! test each sphere to see if it overlaps the sphere of interest ! do i = 1, 2*nres if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30 rplus = rr + radius(i) tx = c(1,i) - xr if (abs(tx) .ge. rplus) goto 30 ty = c(2,i) - yr if (abs(ty) .ge. rplus) goto 30 tz = c(3,i) - zr if (abs(tz) .ge. rplus) goto 30 ! ! check for sphere overlap by testing distance against radii ! xysq = tx*tx + ty*ty if (xysq .lt. delta2) then tx = delta ty = 0.0d0 xysq = delta2 end if ccsq = xysq + tz*tz cc = sqrt(ccsq) if (rplus-cc .le. delta) goto 30 rminus = rr - radius(i) ! ! check to see if sphere of interest is completely buried ! if (cc-abs(rminus) .le. delta) then if (rminus .le. 0.0d0) goto 170 goto 30 end if ! ! check for too many overlaps with sphere of interest ! if (io .ge. maxarc) then write (iout,20) 20 format (/,' SURFATOM -- Increase the Value of MAXARC') stop end if ! ! get overlap between current sphere and sphere of interest ! io = io + 1 xc1(io) = tx yc1(io) = ty zc1(io) = tz dsq1(io) = xysq bsq1(io) = ccsq b1(io) = cc gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io)) intag1(io) = i omit(io) = .false. 30 continue end do ! ! case where no other spheres overlap the sphere of interest ! if (io .eq. 0) then area = 4.0d0 * pi * rrsq return end if ! ! case where only one sphere overlaps the sphere of interest ! if (io .eq. 1) then area = pix2 * (1.0d0 + gr(1)) area = mod(area,4.0d0*pi) * rrsq return end if ! ! case where many spheres intersect the sphere of interest; ! sort the intersecting spheres by their degree of overlap ! call sort2 (io,gr,key) do i = 1, io k = key(i) intag(i) = intag1(k) xc(i) = xc1(k) yc(i) = yc1(k) zc(i) = zc1(k) dsq(i) = dsq1(k) b(i) = b1(k) bsq(i) = bsq1(k) end do ! ! get radius of each overlap circle on surface of the sphere ! do i = 1, io gi = gr(i) * rr bg(i) = b(i) * gi risq(i) = rrsq - gi*gi ri(i) = sqrt(risq(i)) ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i)))) end do ! ! find boundary of inaccessible area on sphere of interest ! do k = 1, io-1 if (.not. omit(k)) then txk = xc(k) tyk = yc(k) tzk = zc(k) bk = b(k) therk = ther(k) ! ! check to see if J circle is intersecting K circle; ! get distance between circle centers and sum of radii ! do j = k+1, io if (omit(j)) goto 60 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j)) cc = acos(min(1.0d0,max(-1.0d0,cc))) td = therk + ther(j) ! ! check to see if circles enclose separate regions ! if (cc .ge. td) goto 60 ! ! check for circle J completely inside circle K ! if (cc+ther(j) .lt. therk) goto 40 ! ! check for circles that are essentially parallel ! if (cc .gt. delta) goto 50 40 continue omit(j) = .true. goto 60 ! ! check to see if sphere of interest is completely buried ! 50 continue if (pix2-cc .le. td) goto 170 60 continue end do end if end do ! ! find T value of circle intersections ! do k = 1, io if (omit(k)) goto 110 omit(k) = .true. narc = 0 top = .false. txk = xc(k) tyk = yc(k) tzk = zc(k) dk = sqrt(dsq(k)) bsqk = bsq(k) bk = b(k) gk = gr(k) * rr risqk = risq(k) rik = ri(k) therk = ther(k) ! ! rotation matrix elements ! t1 = tzk / (bk*dk) axx = txk * t1 axy = tyk * t1 axz = dk / bk ayx = tyk / dk ayy = txk / dk azx = txk / bk azy = tyk / bk azz = tzk / bk do j = 1, io if (.not. omit(j)) then txj = xc(j) tyj = yc(j) tzj = zc(j) ! ! rotate spheres so K vector colinear with z-axis ! uxj = txj*axx + tyj*axy - tzj*axz uyj = tyj*ayy - txj*ayx uzj = txj*azx + tyj*azy + tzj*azz cosine = min(1.0d0,max(-1.0d0,uzj/b(j))) if (acos(cosine) .lt. therk+ther(j)) then dsqj = uxj*uxj + uyj*uyj tb = uzj*gk - bg(j) txb = uxj * tb tyb = uyj * tb td = rik * dsqj tr2 = risqk*dsqj - tb*tb tr2 = max(eps,tr2) tr = sqrt(tr2) txr = uxj * tr tyr = uyj * tr ! ! get T values of intersection for K circle ! tb = (txb+tyr) / td tb = min(1.0d0,max(-1.0d0,tb)) tk1 = acos(tb) if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1 tb = (txb-tyr) / td tb = min(1.0d0,max(-1.0d0,tb)) tk2 = acos(tb) if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j)) if (abs(thec) .lt. 1.0d0) then the = -acos(thec) else if (thec .ge. 1.0d0) then the = 0.0d0 else if (thec .le. -1.0d0) then the = -pi end if ! ! see if "tk1" is entry or exit point; check t=0 point; ! "ti" is exit point, "tf" is entry point ! cosine = min(1.0d0,max(-1.0d0, & (uzj*gk-uxj*rik)/(b(j)*rr))) if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then ti = tk2 tf = tk1 else ti = tk2 tf = tk1 end if narc = narc + 1 if (narc .ge. maxarc) then write (iout,70) 70 format (/,' SURFATOM -- Increase the Value',& ' of MAXARC') stop end if if (tf .le. ti) then arcf(narc) = tf arci(narc) = 0.0d0 tf = pix2 lt(narc) = j ex(narc) = the top = .true. narc = narc + 1 end if arcf(narc) = tf arci(narc) = ti lt(narc) = j ex(narc) = the ux(j) = uxj uy(j) = uyj uz(j) = uzj end if end if end do omit(k) = .false. ! ! special case; K circle without intersections ! if (narc .le. 0) goto 90 ! ! general case; sum up arclength and set connectivity code ! call sort2 (narc,arci,key) arcsum = arci(1) mi = key(1) t = arcf(mi) ni = mi if (narc .gt. 1) then do j = 2, narc m = key(j) if (t .lt. arci(j)) then arcsum = arcsum + arci(j) - t exang = exang + ex(ni) jb = jb + 1 if (jb .ge. maxarc) then write (iout,80) 80 format (/,' SURFATOM -- Increase the Value',& ' of MAXARC') stop end if i = lt(ni) kent(jb) = maxarc*i + k i = lt(m) kout(jb) = maxarc*k + i end if tt = arcf(m) if (tt .ge. t) then t = tt ni = m end if end do end if arcsum = arcsum + pix2 - t if (.not. top) then exang = exang + ex(ni) jb = jb + 1 i = lt(ni) kent(jb) = maxarc*i + k i = lt(mi) kout(jb) = maxarc*k + i end if goto 100 90 continue arcsum = pix2 ib = ib + 1 100 continue arclen = arclen + gr(k)*arcsum 110 continue end do if (arclen .eq. 0.0d0) goto 170 if (jb .eq. 0) goto 150 ! ! find number of independent boundaries and check connectivity ! j = 0 do k = 1, jb if (kout(k) .ne. 0) then i = k 120 continue m = kout(i) kout(i) = 0 j = j + 1 do ii = 1, jb if (m .eq. kent(ii)) then if (ii .eq. k) then ib = ib + 1 if (j .eq. jb) goto 150 goto 130 end if i = ii goto 120 end if end do 130 continue end if end do ib = ib + 1 ! ! attempt to fix connectivity error by moving atom slightly ! if (moved) then write (iout,140) ir 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6) else moved = .true. xr = xr + rmove yr = yr + rmove zr = zr + rmove goto 10 end if ! ! compute the exposed surface area for the sphere of interest ! 150 continue area = ib*pix2 + exang + arclen area = mod(area,4.0d0*pi) * rrsq ! ! attempt to fix negative area by moving atom slightly ! if (area .lt. 0.0d0) then if (moved) then write (iout,160) ir 160 format (/,' SURFATOM -- Negative Area at Atom',i6) else moved = .true. xr = xr + rmove yr = yr + rmove zr = zr + rmove goto 10 end if end if 170 continue return end subroutine surfatom !---------------------------------------------------------------- !---------------------------------------------------------------- subroutine alloc_MD_arrays !EL Allocation of arrays used by MD module integer :: nres2,nres6 nres2=nres*2 nres6=nres*6 !---------------------- #ifndef LANG0 ! commom.langevin ! common /langforc/ allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2) allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6) if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) allocate(fricvec(nres2,nres2)) allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2)) allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2)) allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2) allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch)) allocate(afric0_mat(nres2,nres2,0:maxflag_stoch)) allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch)) allocate(prand0_mat(nres2,nres2,0:maxflag_stoch)) allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch)) allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch) allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch) ! common /langmat/ allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2) !---------------------- #else ! commom.langevin.lang0 ! common /langforc/ allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2) #ifndef FIVEDIAG if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2) #endif allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6) allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch) #endif #ifndef FIVEDIAG if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) #endif !---------------------- ! commom.hairpin in CSA module !---------------------- ! common.mce in MCM_MD module !---------------------- ! common.MD ! common /mdgrad/ in module.energy ! common /back_constr/ in module.energy ! common /qmeas/ in module.energy ! common /mdpar/ ! common /MDcalc/ allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4) ! common /lagrange/ allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2) allocate(d_a_work(nres6)) !(6*MAXRES) #ifdef FIVEDIAG allocate(DM(nres2),DU1(nres2),DU2(nres2)) allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2)) !#ifdef DEBUG allocate(Gvec(1300,1300))!maximum number of protein in one chain !#endif #else write (iout,*) "Before A Allocation",nfgtasks-1 call flush(iout) allocate(Gmat(nres2,nres2),A(nres2,nres2)) if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2) #endif allocate(Geigen(nres2)) !(maxres2) if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2) ! common /inertia/ in io_conf: parmread ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1) ! common /langevin/in io read_MDpar ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1) ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp) ! in io_conf: parmread ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1) ! common /mdpmpi/ in control: ergastulum if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1)) if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1)) if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1) if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs) !---------------------- ! common.muca in read_muca ! common /double_muca/ ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,& ! nemuca2,hist !(4*maxres) ! common /integer_muca/ ! integer :: nmuca,imtime,muca_smooth ! common /mucarem/ ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs) !---------------------- ! common.MD ! common /mdgrad/ in module.energy ! common /back_constr/ in module.energy ! common /qmeas/ in module.energy ! common /mdpar/ ! common /MDcalc/ ! common /lagrange/ allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6)) allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6) allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2) allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES) !---------------------- ! COMMON /BANII/ D allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres ! common /stochcalc/ stochforcvec allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres !---------------------- return end subroutine alloc_MD_arrays !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module MDyn