source/unres/src_MD/src/md-diff/mts/init_md.f

   1 c---------------------------------------------------------
   2       subroutine init_MD
   3 c  Set up the initial conditions of a MD simulation
   4       implicit real*8 (a-h,o-z)
   5       include 'DIMENSIONS'
   6 #ifdef MP
   7       include 'mpif.h'
   8       character*16 form
   9       integer IERROR,ERRCODE
  10 #endif
  11       include 'COMMON.SETUP'
  12       include 'COMMON.CONTROL'
  13       include 'COMMON.VAR'
  14       include 'COMMON.MD'
  15 #ifndef LANG0
  16       include 'COMMON.LANGEVIN'
  17 #else
  18       include 'COMMON.LANGEVIN.lang0'
  19 #endif
  20       include 'COMMON.CHAIN'
  21       include 'COMMON.DERIV'
  22       include 'COMMON.GEO'
  23       include 'COMMON.LOCAL'
  24       include 'COMMON.INTERACT'
  25       include 'COMMON.IOUNITS'
  26       include 'COMMON.NAMES'
  27       include 'COMMON.REMD'
  28       real*8 energia_long(0:n_ene),
  29      &  energia_short(0:n_ene),vcm(3),incr(3)
  30       double precision cm(3),L(3),xv,sigv,lowb,highb
  31       double precision varia(maxvar)
  32       character*256 qstr
  33       integer ilen
  34       external ilen
  35       character*50 tytul
  36       logical file_exist
  37       common /gucio/ cm
  38       d_time0=d_time
  39 c      write(iout,*) "d_time", d_time
  40 c Compute the standard deviations of stochastic forces for Langevin dynamics
  41 c if the friction coefficients do not depend on surface area
  42       if (lang.gt.0 .and. .not.surfarea) then
  43         do i=nnt,nct-1
  44           stdforcp(i)=stdfp*dsqrt(gamp)
  45         enddo
  46         do i=nnt,nct
  47           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
  48         enddo
  49       endif
  50 c Open the pdb file for snapshotshots
  51 #ifdef MPI
  52       if(mdpdb) then
  53         if (ilen(tmpdir).gt.0)
  54      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
  55      &      liczba(:ilen(liczba))//".pdb")
  56         open(ipdb,
  57      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
  58      &  //".pdb")
  59       else
  60 #ifdef NOXDR
  61         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
  62      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
  63      &      liczba(:ilen(liczba))//".x")
  64         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
  65      &  //".x"
  66 #else
  67         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
  68      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
  69      &      liczba(:ilen(liczba))//".cx")
  70         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
  71      &  //".cx"
  72 #endif
  73       endif
  74 #else
  75       if(mdpdb) then
  76          if (ilen(tmpdir).gt.0)
  77      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
  78          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
  79       else
  80          if (ilen(tmpdir).gt.0)
  81      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
  82          cartname=prefix(:ilen(prefix))//"_MD.cx"
  83       endif
  84 #endif
  85       if (usampl) then
  86         write (qstr,'(256(1h ))')
  87         ipos=1
  88         do i=1,nfrag
  89           iq = qinfrag(i,iset)*10
  90           iw = wfrag(i,iset)/100
  91           if (iw.gt.0) then
  92             if(me.eq.king.or..not.out1file)
  93      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
  94             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
  95             ipos=ipos+7
  96           endif
  97         enddo
  98         do i=1,npair
  99           iq = qinpair(i,iset)*10
 100           iw = wpair(i,iset)/100
 101           if (iw.gt.0) then
 102             if(me.eq.king.or..not.out1file)
 103      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
 104             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
 105             ipos=ipos+7
 106           endif
 107         enddo
 108 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
 109 #ifdef NOXDR
 110 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
 111 #else
 112 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
 113 #endif
 114 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
 115       endif
 116       icg=1
 117       if (rest) then
 118        if (restart1file) then
 119          if (me.eq.king)
 120      &     inquire(file=mremd_rst_name,exist=file_exist)
 121            write (*,*) me," Before broadcast: file_exist",file_exist
 122          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
 123      &          IERR)
 124          write (*,*) me," After broadcast: file_exist",file_exist
 125 c        inquire(file=mremd_rst_name,exist=file_exist)
 126         if(me.eq.king.or..not.out1file)
 127      &   write(iout,*) "Initial state read by master and distributed"
 128        else
 129          if (ilen(tmpdir).gt.0)
 130      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
 131      &      //liczba(:ilen(liczba))//'.rst')
 132         inquire(file=rest2name,exist=file_exist)
 133        endif
 134        if(file_exist) then
 135          if(.not.restart1file) then
 136            if(me.eq.king.or..not.out1file)
 137      &      write(iout,*) "Initial state will be read from file ",
 138      &      rest2name(:ilen(rest2name))
 139            call readrst
 140          endif
 141          call rescale_weights(t_bath)
 142        else
 143         if(me.eq.king.or..not.out1file)then
 144          if (restart1file) then
 145           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
 146      &       " does not exist"
 147          else
 148           write(iout,*) "File ",rest2name(:ilen(rest2name)),
 149      &       " does not exist"
 150          endif
 151          write(iout,*) "Initial velocities randomly generated"
 152         endif
 153         call random_vel
 154         totT=0.0d0
 155        endif
 156       else
 157 c Generate initial velocities
 158         if(me.eq.king.or..not.out1file)
 159      &   write(iout,*) "Initial velocities randomly generated"
 160         call random_vel
 161         totT=0.0d0
 162       endif
 163 c      rest2name = prefix(:ilen(prefix))//'.rst'
 164       if(me.eq.king.or..not.out1file)then
 165        write (iout,*) "Initial velocities"
 166        do i=0,nres
 167          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 168      &   (d_t(j,i+nres),j=1,3)
 169        enddo
 170 c  Zeroing the total angular momentum of the system
 171        write(iout,*) "Calling the zero-angular
 172      & momentum subroutine"
 173       endif
 174       call inertia_tensor
 175 c  Getting the potential energy and forces and velocities and accelerations
 176       call vcm_vel(vcm)
 177 c      write (iout,*) "velocity of the center of the mass:"
 178 c      write (iout,*) (vcm(j),j=1,3)
 179       do j=1,3
 180         d_t(j,0)=d_t(j,0)-vcm(j)
 181       enddo
 182 c Removing the velocity of the center of mass
 183       call vcm_vel(vcm)
 184       if(me.eq.king.or..not.out1file)then
 185        write (iout,*) "vcm right after adjustment:"
 186        write (iout,*) (vcm(j),j=1,3)
 187       endif
 188       if (.not.rest) then
 189          call chainbuild
 190          if(iranconf.ne.0) then
 191           if (overlapsc) then
 192            print *, 'Calling OVERLAP_SC'
 193            call overlap_sc(fail)
 194           endif
 195
 196           if (searchsc) then
 197            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
 198            print *,'SC_move',nft_sc,etot
 199            if(me.eq.king.or..not.out1file)
 200      &      write(iout,*) 'SC_move',nft_sc,etot
 201           endif
 202
 203           if(dccart)then
 204            print *, 'Calling MINIM_DC'
 205            call minim_dc(etot,iretcode,nfun)
 206           else
 207            call geom_to_var(nvar,varia)
 208            print *,'Calling MINIMIZE.'
 209            call minimize(etot,varia,iretcode,nfun)
 210            call var_to_geom(nvar,varia)
 211           endif
 212           if(me.eq.king.or..not.out1file)
 213      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
 214          endif
 215       endif
 216       call chainbuild_cart
 217       call kinetic(EK)
 218       if (tbf) then
 219         call verlet_bath(EK)
 220       endif
 221       kinetic_T=2.0d0/(dimen3*Rb)*EK
 222       if(me.eq.king.or..not.out1file)then
 223        call cartprint
 224        call intout
 225       endif
 226 #ifdef MPI
 227       tt0=MPI_Wtime()
 228 #else
 229       tt0=tcpu()
 230 #endif
 231       call zerograd
 232       call etotal(potEcomp)
 233 #ifdef TIMING_ENE
 234 #ifdef MPI
 235       t_etotal=t_etotal+MPI_Wtime()-tt0
 236 #else
 237       t_etotal=t_etotal+tcpu()-tt0
 238 #endif
 239 #endif
 240       potE=potEcomp(0)
 241       call cartgrad
 242       call lagrangian
 243       call max_accel
 244       if (amax*d_time .gt. dvmax) then
 245         d_time=d_time*dvmax/amax
 246         if(me.eq.king.or..not.out1file) write (iout,*)
 247      &   "Time step reduced to",d_time,
 248      &   " because of too large initial acceleration."
 249       endif
 250       if(me.eq.king.or..not.out1file)then
 251        write(iout,*) "Potential energy and its components"
 252        call enerprint(potEcomp)
 253 c       write(iout,*) (potEcomp(i),i=0,n_ene)
 254       endif
 255       potE=potEcomp(0)-potEcomp(20)
 256       totE=EK+potE
 257       itime=0
 258       if (ntwe.ne.0) call statout(itime)
 259       if(me.eq.king.or..not.out1file)
 260      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
 261      &   " Kinetic energy",EK," potential energy",potE,
 262      &   " total energy",totE," maximum acceleration ",
 263      &   amax
 264       if (large) then
 265         write (iout,*) "Initial coordinates"
 266         do i=1,nres
 267           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
 268      &    (c(j,i+nres),j=1,3)
 269         enddo
 270         write (iout,*) "Initial dC"
 271         do i=0,nres
 272           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
 273      &    (dc(j,i+nres),j=1,3)
 274         enddo
 275         write (iout,*) "Initial velocities"
 276         write (iout,"(13x,' backbone ',23x,' side chain')")
 277         do i=0,nres
 278           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 279      &    (d_t(j,i+nres),j=1,3)
 280         enddo
 281         write (iout,*) "Initial accelerations"
 282         do i=0,nres
 283 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 284           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
 285      &    (d_a(j,i+nres),j=1,3)
 286         enddo
 287       endif
 288       do i=0,2*nres
 289         do j=1,3
 290           dc_old(j,i)=dc(j,i)
 291           d_t_old(j,i)=d_t(j,i)
 292           d_a_old(j,i)=d_a(j,i)
 293         enddo
 294 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
 295       enddo
 296       if (RESPA) then
 297 #ifdef MPI
 298       tt0 =MPI_Wtime()
 299 #else
 300       tt0 = tcpu()
 301 #endif
 302         call zerograd
 303         call etotal_short(energia_short)
 304 #ifdef TIMING_ENE
 305 #ifdef MPI
 306         t_eshort=t_eshort+MPI_Wtime()-tt0
 307 #else
 308         t_eshort=t_eshort+tcpu()-tt0
 309 #endif
 310 #endif
 311         call cartgrad
 312         call lagrangian
 313         if(.not.out1file .and. large) then
 314           write (iout,*) "energia_long",energia_long(0),
 315      &     " energia_short",energia_short(0),
 316      &     " total",energia_long(0)+energia_short(0)
 317           write (iout,*) "Initial fast-force accelerations"
 318           do i=0,nres
 319             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 320      &      (d_a(j,i+nres),j=1,3)
 321           enddo
 322         endif
 323 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
 324         do i=0,2*nres
 325           do j=1,3
 326             d_a_short(j,i)=d_a(j,i)
 327           enddo
 328         enddo
 329 #ifdef MPI
 330         tt0=MPI_Wtime()
 331 #else
 332         tt0=tcpu()
 333 #endif
 334         call zerograd
 335         call etotal_long(energia_long)
 336 #ifdef TIMING_ENE
 337 #ifdef MPI
 338         t_elong=t_elong+MPI_Wtime()-tt0
 339 #else
 340         t_elong=t_elong+tcpu()-tt0
 341 #endif
 342 #endif
 343         call cartgrad
 344         call lagrangian
 345         if(.not.out1file .and. large) then
 346           write (iout,*) "energia_long",energia_long(0)
 347           write (iout,*) "Initial slow-force accelerations"
 348           do i=0,nres
 349             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 350      &      (d_a(j,i+nres),j=1,3)
 351           enddo
 352         endif
 353 #ifdef MPI
 354         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 355 #else
 356         t_enegrad=t_enegrad+tcpu()-tt0
 357 #endif
 358       endif
 359       return
 360       end