X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fsrc%2Fmd-diff%2Fmts%2Finit_md.f;fp=source%2Funres%2Fsrc_MD%2Fsrc%2Fmd-diff%2Fmts%2Finit_md.f;h=0000000000000000000000000000000000000000;hb=0a11a2c4ccee14ed99ae44f2565b270ba8d4bbb6;hp=05e2460092c6e49293a8e139771895dfc767a728;hpb=5eb407964903815242c59de10960f42761139e10;p=unres.git diff --git a/source/unres/src_MD/src/md-diff/mts/init_md.f b/source/unres/src_MD/src/md-diff/mts/init_md.f deleted file mode 100644 index 05e2460..0000000 --- a/source/unres/src_MD/src/md-diff/mts/init_md.f +++ /dev/null @@ -1,360 +0,0 @@ -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