+++ /dev/null
-c---------------------------------------------------------
- subroutine init_MD
-c Set up the initial conditions of a MD simulation
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
-#ifdef MP
- include 'mpif.h'
- include 'COMMON.INFO'
- include 'COMMON.SETUP'
- character*4 liczba
- character*16 form
-#endif
- include 'COMMON.CONTROL'
- include 'COMMON.VAR'
- include 'COMMON.MD'
- include 'COMMON.LANGEVIN'
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.NAMES'
- real*8 energia(0:n_ene),energia_long(0:n_ene),
- & energia_short(0:n_ene),vcm(3),incr(3),E_short
- double precision cm(3),L(3),xv,sigv,lowb,highb
- character*256 qstr
- integer ilen
- external ilen
- character*50 tytul
- common /gucio/ cm
- d_time0=d_time
- write(iout,*) "d_time", d_time
-c Compute the standard deviations of stochastic forces for Langevin dynamics
-c 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(itype(i))*dsqrt(gamsc(itype(i)))
- enddo
- endif
-c Open the pdb file for snapshotshots
-#ifdef MPI
- if (nprocs.eq.1) then
- npos=3
- else
- npos = dlog10(dfloat(nprocs-1))+1
- endif
- if (npos.lt.3) npos=3
- write (liczba,'(i1)') npos
- form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
- & //')'
- write (liczba,form) myrank
- if(mdpdb) then
- open(ipdb,
- & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
- & //".pdb")
- else
- cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
- & //".cx"
- endif
-#else
- if(mdpdb) then
- open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
- else
- cartname=prefix(:ilen(prefix))//"_MD.cx"
- endif
-#endif
- if (usampl) then
- write (qstr,'(256(1h ))')
- ipos=1
- do i=1,nfrag
- iq = qinfrag(i)*10
- iw = wfrag(i)/100
- if (iw.gt.0) then
- write (iout,*) "Frag",qinfrag(i),wfrag(i),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)*10
- iw = wpair(i)/100
- if (iw.gt.0) then
- write (iout,*) "Pair",i,qinpair(i),wpair(i),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'
- cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
- statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
- endif
- icg=1
- if (rest) then
- write(iout,*) "Initial state will be read from file ",
- & rest2name(:ilen(rest2name))
- call readrst
- else
-c Generate initial velocities
- write(iout,*) "Initial velocities randomly generated"
- call random_vel
- totT=0.0d0
- endif
-c rest2name = prefix(:ilen(prefix))//'.rst'
- write(iout,*) "Initial backbone velocities"
- do i=nnt,nct-1
- write(iout,*) (d_t(j,i),j=1,3)
- enddo
- write(iout,*) "Initial side-chain velocities"
- do i=nnt,nct
- write(iout,*) (d_t(j,i+nres),j=1,3)
- enddo
-c Zeroing the total angular momentum of the system
- write(iout,*) "Calling the zero-angular
- & momentum subroutine"
- call inertia_tensor
-c Getting the potential energy and forces and velocities and accelerations
- call vcm_vel(vcm)
-c write (iout,*) "velocity of the center of the mass:"
-c write (iout,*) (vcm(j),j=1,3)
- do j=1,3
- d_t(j,0)=d_t(j,0)-vcm(j)
- enddo
-c Removing the velocity of the center of mass
- call vcm_vel(vcm)
- write (iout,*) "vcm right after adjustment:"
- write (iout,*) (vcm(j),j=1,3)
- if (.not.rest) then
- call chainbuild
- endif
- call chainbuild_cart
- call kinetic(EK)
- if (tbf) then
- call verlet_bath(EK)
- endif
- kinetic_T=2.0d0/(dimen*Rb)*EK
- call cartprint
- call intout
- call zerograd
- call etotal(energia)
- potE=energia(0)
-
- if(tnp .or. tnp1) then
- s_np=1.0
- pi_np=0.0
- HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
- H0=Hnose1
- write(iout,*) 'H0= ',H0
- endif
-
- if(tnh) then
- HNose1=Hnose_nh(EK,potE)
- H0=HNose1
- write (iout,*) 'H0= ',H0
- endif
-
-
- call cartgrad
- call lagrangian
- call max_accel
- if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
- write(iout,*) "Potential energy"
- write(iout,*) (energia(i),i=0,n_ene)
- potE=energia(0)-energia(20)
- totE=EK+potE
- itime=0
- call statout(itime)
- write (iout,*) "Initial:",
- & " Kinetic energy",EK," potential energy",potE,
- & " total energy",totE," maximum acceleration ",
- & amax
- if (large) 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
- write (iout,*) "Initial 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
- 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
-c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
- enddo
- if (RESPA) then
- call zerograd
- call etotal_long(energia_long)
- E_long=energia_long(0)
- if(tnp .or. tnp1) then
- call etotal_short(energia_short)
- E_short=energia_short(0)
- HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen)
- Csplit=Hnose1
-c Csplit =110
-c_new_var_csplit Csplit=H0-E_long
-c Csplit = H0-energia_short(0)
- write(iout,*) 'Csplit= ',Csplit
- endif
- call cartgrad
- call lagrangian
-c call etotal_short(energia_short)
-c write (iout,*) "energia_long",energia_long(0),
-c & " energia_short",energia_short(0),
-c & " total",energia_long(0)+energia_short(0)
- endif
- return
- end
projects / unres.git / blobdiff