+++ /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'
- character*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*8 energia_long(0:n_ene),
- & energia_short(0:n_ene),vcm(3),incr(3)
- double precision cm(3),L(3),xv,sigv,lowb,highb
- double precision varia(maxvar)
- character*256 qstr
- integer ilen
- external ilen
- character*50 tytul
- logical file_exist
- common /gucio/ cm
- d_time0=d_time
-c 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(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
-c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
-#ifdef NOXDR
-c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
-#else
-c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
-#endif
-c 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
- call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
- & IERR)
- write (*,*) me," After broadcast: file_exist",file_exist
-c 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
-c Generate initial velocities
- if(me.eq.king.or..not.out1file)
- & write(iout,*) "Initial velocities randomly generated"
- call random_vel
- totT=0.0d0
- endif
-c 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
-c Zeroing the total angular momentum of the system
- write(iout,*) "Calling the zero-angular
- & momentum subroutine"
- endif
- 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)
- 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(EK)
- 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)
-#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)
-c 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
-c 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
-c 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)
-#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
-C 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)
-#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
projects / unres.git / blobdiff