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 !el real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres) !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- 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(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(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres nbond=nct-nnt do i=nnt,nct if (itype(i).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).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).ne.10) then ind=ind+1 Td(i)=Td(i)+vbldsc0(1,itype(k))*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).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).ne.10) then ind=ind+1 xx=vbld(i+nres)-vbldsc0(1,itype(i)) 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).ne.10) then ind=ind+1 blen2 = scalar(dc(1,i+nres),dc(1,i+nres)) ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2 diffbond=dabs(vbldsc0(1,itype(i))-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).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).ne.10) then ind=ind+1 xx=vbld(i+nres)-vbldsc0(1,itype(i)) 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(20) call cartgrad totT=totT+d_time ! 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 !----------------------------------------------------------------------------- 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 integer :: i,j,k,iti 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)),i=1,nres) ! The translational part for peptide virtual bonds do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct-1 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3) do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) enddo vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) KEt_p=KEt_p+(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 iti=iabs(itype(i)) if (itype(i).eq.10) 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)," v",(v(j),j=1,3) KEt_sc=KEt_sc+msc(iti)*(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 ! 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+(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 iti=iabs(itype(i)) if (itype(i).ne.10) 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)*(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*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc) ! write (iout,*) "KE_total",KE_total return end subroutine kinetic !----------------------------------------------------------------------------- ! 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 #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 :: itime,i,j,nharp integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3) ! logical :: ovrtim real(kind=8) :: tt0,scalfac integer :: nres2 nres2=2*nres ! #ifdef MPI 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 ! Determine the inverse of the inertia matrix. call setup_MD_matrices ! Initialize MD call 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 call 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 (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 call 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 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. call velverlet_step(itime) 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 if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0) call statout(itime) #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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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 write (tytul,'("time",f8.2)') totT if(mdpdb) then call hairpin(.true.,nharp,iharp) call secondary2(.true.) call pdbout(potE,tytul,ipdb) else call cartout(totT) 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 ! 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 #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 = 20 !el common /gucio/ cm !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec integer :: itime,icount_scale,itime_scal,i,j,ifac_time logical :: scale real(kind=8) :: epdrift,tt0,fac_time ! if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres scale=.true. icount_scale=0 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 itime_scal=0 do while (scale) icount_scale=icount_scale+1 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.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) 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(20) 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) if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step scale=.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 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 scale=.false. 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 :: itime,itt,i,j,itsplit 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) 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 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 (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 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(20) ! potE=energia_short(0)+energia_long(0) totT=totT+d_time ! 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 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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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 #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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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 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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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 ! ! Compute friction and stochastic forces ! #ifdef MPI time00=MPI_Wtime() call 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) 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 ! ! 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 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 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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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 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(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0 integer :: i,j,ind,inres ! 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) ! call ginv_mult(stochforcvec, d_as_work1) ! ! 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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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 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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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).ne.10 .and. itype(i).ne.ntyp1) 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 ! 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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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 ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MP include 'mpif.h' character(len=16) :: form 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.REMD' real(kind=8),dimension(0:n_ene) :: energia_long,energia_short 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,nfun,itime,ierr 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 stdforcp(i)=stdfp*dsqrt(gamp) enddo do i=nnt,nct stdforcsc(i)=stdfsc(iabs(itype(i))) & *dsqrt(gamsc(iabs(itype(i)))) 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 endif else ! Generate initial velocities if(me.eq.king.or..not.out1file) & write(iout,*) "Initial velocities randomly generated" call random_vel totT=0.0d0 endif ! rest2name = prefix(:ilen(prefix))//'.rst' if(me.eq.king.or..not.out1file)then write (iout,*) "Initial 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 ! 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 if (.not.rest) then call chainbuild if(iranconf.ne.0) then if (overlapsc) then print *, 'Calling OVERLAP_SC' call overlap_sc(fail) endif 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) 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 endif endif call chainbuild_cart call kinetic(EK) if (tbf) then call verlet_bath endif kinetic_T=2.0d0/(dimen3*Rb)*EK if(me.eq.king.or..not.out1file)then call cartprint call intout endif #ifdef MPI tt0=MPI_Wtime() #else tt0=tcpu() #endif call zerograd 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 call lagrangian 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(20) totE=EK+potE itime=0 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,'(i3,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 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs real(kind=8) :: sumx #ifdef DEBUG real(kind=8) ,allocatable, dimension(:) :: rsold real (kind=8),allocatable,dimension(:,:) :: matold #endif #endif integer :: i,j,ii,k,ind,mark,imark ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m ! First generate velocities in the eigenspace of the G matrix ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3 ! call flush(iout) 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 #ifdef FIVEDIAG ! Transform velocities to UNRES coordinate space allocate (DL1(2*nres)) allocate (DDU1(2*nres)) allocate (DL2(2*nres)) allocate (DDU2(2*nres)) allocate (xsolv(2*nres)) allocate (DML(2*nres)) allocate (rs(2*nres)) #ifdef DEBUG allocate (rsold(2*nres)) allocate (matold(2*nres,2*nres)) matold=0 matold(1,1)=DMorig(1) matold(1,2)=DU1orig(1) matold(1,3)=DU2orig(1) write (*,*) DMorig(1),DU1orig(1),DU2orig(1) do i=2,dimen-1 do j=1,dimen if (i.eq.j) then matold(i,j)=DMorig(i) matold(i,j-1)=DU1orig(i-1) if (j.ge.3) then matold(i,j-2)=DU2orig(i-2) endif endif enddo do j=1,dimen-1 if (i.eq.j) then matold(i,j+1)=DU1orig(i) end if enddo do j=1,dimen-2 if(i.eq.j) then matold(i,j+2)=DU2orig(i) endif enddo enddo matold(dimen,dimen)=DMorig(dimen) matold(dimen,dimen-1)=DU1orig(dimen-1) matold(dimen,dimen-2)=DU2orig(dimen-2) write(iout,*) "old gmatrix" call matout(dimen,dimen,2*nres,2*nres,matold) #endif d_t_work=0.0d0 do i=1,dimen ! Find the ith eigenvector of the pentadiagonal inertiq matrix do imark=dimen,1,-1 do j=1,imark-1 DML(j)=DMorig(j)-geigen(i) enddo do j=imark+1,dimen DML(j-1)=DMorig(j)-geigen(i) enddo do j=1,imark-2 DDU1(j)=DU1orig(j) enddo DDU1(imark-1)=DU2orig(imark-1) do j=imark+1,dimen-1 DDU1(j-1)=DU1orig(j) enddo do j=1,imark-3 DDU2(j)=DU2orig(j) enddo DDU2(imark-2)=0.0d0 DDU2(imark-1)=0.0d0 do j=imark,dimen-3 DDU2(j)=DU2orig(j+1) enddo do j=1,dimen-3 DL2(j+2)=DDU2(j) enddo do j=1,dimen-2 DL1(j+1)=DDU1(j) enddo #ifdef DEBUG write (iout,*) "DL2,DL1,DML,DDU1,DDU2" write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1) write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1) write(iout,'(10f10.5)') (DML(k),k=1,dimen-1) write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2) write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3) #endif rs=0.0d0 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2) if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1) if (imark.lt.dimen) rs(imark)=-DU1orig(imark) if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark) #ifdef DEBUG rsold=rs #endif ! write (iout,*) "Vector rs" ! do j=1,dimen-1 ! write (iout,*) j,rs(j) ! enddo xsolv=0.0d0 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark) if (mark.eq.1) then #ifdef DEBUG ! Check solution do j=1,imark-1 sumx=-geigen(i)*xsolv(j) do k=1,imark-1 sumx=sumx+matold(j,k)*xsolv(k) enddo do k=imark+1,dimen sumx=sumx+matold(j,k)*xsolv(k-1) enddo write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j) enddo do j=imark+1,dimen sumx=-geigen(i)*xsolv(j-1) do k=1,imark-1 sumx=sumx+matold(j,k)*xsolv(k) enddo do k=imark+1,dimen sumx=sumx+matold(j,k)*xsolv(k-1) enddo write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1) enddo ! end check write (iout,*)& "Solution of equations system",i," for eigenvalue",geigen(i) do j=1,dimen-1 write(iout,'(i5,f10.5)') j,xsolv(j) enddo #endif do j=dimen-1,imark,-1 xsolv(j+1)=xsolv(j) enddo xsolv(imark)=1.0d0 #ifdef DEBUG write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i) do j=1,dimen write(iout,'(i5,f10.5)') j,xsolv(j) enddo #endif ! Normalize ith eigenvector sumx=0.0d0 do j=1,dimen sumx=sumx+xsolv(j)**2 enddo sumx=dsqrt(sumx) do j=1,dimen xsolv(j)=xsolv(j)/sumx enddo #ifdef DEBUG write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i) do j=1,dimen write(iout,'(i5,3f10.5)') j,xsolv(j) enddo #endif ! All done at this point for eigenvector i, exit loop exit endif ! mark.eq.1 enddo ! imark if (mark.ne.1) then write (iout,*) "Unable to find eigenvector",i endif ! write (iout,*) "i=",i do k=1,3 ! write (iout,*) "k=",k do j=1,dimen ind=(j-1)*3+k ! write(iout,*) "ind",ind," ind1",3*(i-1)+k d_t_work(ind)=d_t_work(ind) & +xsolv(j)*d_t_work_new(3*(i-1)+k) enddo enddo enddo ! i (loop over the eigenvectors) #ifdef DEBUG write (iout,*) "d_t_work" do i=1,3*dimen write (iout,"(i5,f10.5)") i,d_t_work(i) enddo Ek1=0.0d0 ii=0 do i=nnt,nct if (itype(i).eq.10) then j=ii+3 else j=ii+6 endif if (i.lt.nct) then do k=1,3 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1 Ek1=Ek1+0.5d0*mp*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+& 0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2 enddo endif if (itype(i).ne.10) ii=ii+3 write (iout,*) "i",i," itype",itype(i)," mass",msc(itype(i)) write (iout,*) "ii",ii do k=1,3 ii=ii+1 write (iout,*) "k",k," ii",ii,"EK1",EK1 if (iabs(itype(i)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i)))*(d_t_work(ii)-d_t_work(ii-3))**2 Ek1=Ek1+0.5d0*msc(iabs(itype(i)))*d_t_work(ii)**2 enddo write (iout,*) "i",i," ii",ii enddo write (iout,*) "Ek from d_t_work",Ek1 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) #endif deallocate(DDU1) deallocate(DDU2) deallocate(DL2) deallocate(DL1) deallocate(xsolv) deallocate(DML) deallocate(rs) #ifdef DEBUG deallocate(matold) deallocate(rsold) #endif ind=1 do j=nnt,nct do k=1,3 d_t(k,j)=d_t_work(ind) ind=ind+1 enddo if (itype(j).ne.10 .and. itype(j).ne.ntyp1) then do k=1,3 d_t(k,j+nres)=d_t_work(ind) ind=ind+1 enddo end if enddo #ifdef DEBUG write (iout,*) "Random velocities in the Calpha,SC space" do i=nnt,nct write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3) enddo do i=nnt,nct write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3) enddo #endif do j=1,3 d_t(j,0)=d_t(j,nnt) enddo do i=nnt,nct if (itype(i).eq.10) 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 #ifdef DEBUG write (iout,*)"Random velocities in the virtual-bond-vector space" do i=nnt,nct-1 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3) enddo do i=nnt,nct write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3) enddo call kinetic(Ek1) write (iout,*) "Ek from d_t_work",Ek1 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) #endif #else 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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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) return #undef DEBUG 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 if (itype(i).ne.10 .and. itype(i).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 if (itype(i).ne.10) 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 if (itype(i).ne.10 .and. itype(i).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 !----------------------------------------------------------------------------- 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).ne.10 .and. itype(i).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 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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 if (itype(i).ne.10 .and. itype(i).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 !----------------------------------------------------------------------------- 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 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 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 do i=nnt,nct-1 do j=1,3 cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i) enddo enddo do j=1,3 cm(j)=mp*cm(j) enddo M_SC=0.0d0 do i=nnt,nct iti=iabs(itype(i)) M_SC=M_SC+msc(iabs(iti)) inres=i+nres do j=1,3 cm(j)=cm(j)+msc(iabs(iti))*c(j,inres) enddo enddo do j=1,3 cm(j)=cm(j)/(M_SC+(nct-nnt)*mp) enddo do i=nnt,nct-1 do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3)) Im(1,2)=Im(1,2)-mp*pr(1)*pr(2) Im(1,3)=Im(1,3)-mp*pr(1)*pr(3) Im(2,3)=Im(2,3)-mp*pr(2)*pr(3) Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1)) Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2)) enddo do i=nnt,nct iti=iabs(itype(i)) inres=i+nres do j=1,3 pr(j)=c(j,inres)-cm(j) enddo Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3)) Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2) Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3) Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3) Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1)) Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2)) enddo do i=nnt,nct-1 Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))* & vbld(i+1)*vbld(i+1)*0.25d0 Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 enddo do i=nnt,nct if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then iti=iabs(itype(i)) inres=i+nres Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* & dc_norm(1,inres))*vbld(inres)*vbld(inres) Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)* & dc_norm(2,inres))*vbld(inres)*vbld(inres) Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)* & dc_norm(2,inres))*vbld(inres)*vbld(inres) Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) endif enddo call angmom(cm,L) ! write(iout,*) "The angular momentum before adjustment:" ! write(iout,*) (L(j),j=1,3) 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 d_t(j,i)=d_t(j,i)-vp(j) enddo enddo do i=nnt,nct if(itype(i).ne.10 .and. itype(i).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),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp integer :: iti,inres,i,j ! 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 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*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*vp(j) enddo enddo do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct iti=iabs(itype(i)) inres=i+nres do j=1,3 pr(j)=c(j,inres)-cm(j) enddo if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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) ! 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)+msc(iabs(iti))*vp(j) enddo ! write (iout,*) "L",(l(j),j=1,3) if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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)*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 do j=1,3 vcm(j)=0.0d0 vv(j)=d_t(j,0) enddo summas=0.0d0 do i=nnt,nct if (i.lt.nct) then summas=summas+mp do j=1,3 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i)) enddo endif amas=msc(iabs(itype(i))) summas=summas+amas if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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 !----------------------------------------------------------------------------- ! 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 if (itype(i).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)=Ginv(ind,ind1) enddo enddo do k=nnt,nct if (itype(k).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 enddo do i=nnt,nct if (itype(i).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)=Ginv(ind,ind1) enddo enddo do k=nnt,nct if (itype(k).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).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).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).ne.10) then ind=ind+1 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**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 if (itype(i).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 ! 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).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 if (itype(i).ne.10) then ind=ind+1 xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**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).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 if (itype(k).ne.10) 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 if (itype(i).ne.10) 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 if (itype(i).ne.10) 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 if (itype(i).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),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 if (itype(i).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 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 if (itype(i).ne.10) then 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 = .true., lprn1 = .true., 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).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).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).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).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).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).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).ne.10) then ind=ind+1 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**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).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).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).ne.10) then ind=ind+1 xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**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).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 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 = .false., checkmode = .false. integer :: i,j,ind,k,nres2,nres6 nres2=2*nres nres6=6*nres 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 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) 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 if ((itype(i).ne.10).and.(itype(i).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 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 logical :: lprn = .false. real(kind=8) :: dtdi !el ,gamvec(2*nres) !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy real(kind=8),dimension(2*nres,2*nres) :: 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 #ifdef MPI if (fg_rank.ne.king) goto 10 #endif nres2=2*nres nres6=6*nres if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2) if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres !el 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 !el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) ! 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 ind1=ind1+1 gamvec(ind1)=gamp enddo ! Load the friction coefficients corresponding to side chains m=nct-nnt ind=0 gamsc(ntyp1)=1.0d0 do i=nnt,nct ind=ind+1 ii = ind+m iti=itype(i) gamvec(ii)=gamsc(iabs(iti)) 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 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 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 ! 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).ne.ntyp1) 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 if ((itype(i).ne.10).and.(itype(i).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 if ((itype(i).ne.10).and.(itype(i).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 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 ! ! 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 radius(i)=pstok enddo ! Load side chain radii do i=nnt,nct iti=itype(i) radius(i+nres)=restok(iti) 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 stdforcp(i)=stdfp*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 stdforcsc(i)=stdfsc(itype(i))*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) if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2) allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6) allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch) #endif !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !---------------------- ! 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)) #else 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