source/unres/src_MD-M/MD_A-MTS.F

   1       subroutine MD
   2 c------------------------------------------------
   3 c  The driver for molecular dynamics subroutines
   4 c------------------------------------------------
   5       implicit real*8 (a-h,o-z)
   6       include 'DIMENSIONS'
   7 #ifdef MPI
   8       include "mpif.h"
   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.TIME1'
  28       include 'COMMON.HAIRPIN'
  29       double precision cm(3),L(3),vcm(3)
  30 #ifdef VOUT
  31       double precision v_work(maxres6),v_transf(maxres6)
  32 #endif
  33       integer ilen,rstcount
  34       external ilen
  35       character*50 tytul
  36       common /gucio/ cm
  37       integer itime
  38       logical ovrtim
  39       integer nharp,iharp(4,maxres/3)
  40 c
  41 #ifdef MPI
  42       if (ilen(tmpdir).gt.0)
  43      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
  44      &        //liczba(:ilen(liczba))//'.rst')
  45 #else
  46       if (ilen(tmpdir).gt.0)
  47      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
  48 #endif
  49       t_MDsetup=0.0d0
  50       t_langsetup=0.0d0
  51       t_MD=0.0d0
  52       t_enegrad=0.0d0
  53       t_sdsetup=0.0d0
  54       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
  55       write (iout,'(/a)')
  56      &  "Cartesian coordinates of the initial structure"
  57       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
  58      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
  59       do ires=1,nres
  60         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
  61      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
  62      &    (c(j,ires+nres),j=1,3)
  63       enddo
  64       write (iout,'(/a)')
  65      &  "Initial dC vectors of the chain"
  66       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
  67      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
  68       do ires=1,nres
  69         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
  70      &    restyp(itype(ires)),ires,(dc(j,ires),j=1,3),
  71      &    (dc(j,ires+nres),j=1,3)
  72       enddo
  73       write (iout,'(/a)')
  74      &  "Initial dC_norm vectors of the chain"
  75       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
  76      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
  77       do ires=1,nres
  78         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
  79      &    restyp(itype(ires)),ires,(dc_norm(j,ires),j=1,3),
  80      &    (dc_norm(j,ires+nres),j=1,3)
  81       enddo
  82 #ifdef MPI
  83       tt0=MPI_Wtime()
  84 #else
  85       tt0 = tcpu()
  86 #endif
  87 c Determine the inverse of the inertia matrix.
  88       call setup_MD_matrices
  89 c Initialize MD
  90       call init_MD
  91 #ifdef MPI
  92       t_MDsetup = MPI_Wtime()-tt0
  93 #else
  94       t_MDsetup = tcpu()-tt0
  95 #endif
  96       rstcount=0
  97 c   Entering the MD loop
  98 #ifdef MPI
  99       tt0 = MPI_Wtime()
 100 #else
 101       tt0 = tcpu()
 102 #endif
 103       if (lang.eq.2 .or. lang.eq.3) then
 104 #ifndef   LANG0
 105         call setup_fricmat
 106         if (lang.eq.2) then
 107           call sd_verlet_p_setup
 108         else
 109           call sd_verlet_ciccotti_setup
 110         endif
 111         do i=1,dimen
 112           do j=1,dimen
 113             pfric0_mat(i,j,0)=pfric_mat(i,j)
 114             afric0_mat(i,j,0)=afric_mat(i,j)
 115             vfric0_mat(i,j,0)=vfric_mat(i,j)
 116             prand0_mat(i,j,0)=prand_mat(i,j)
 117             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 118             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 119           enddo
 120         enddo
 121         flag_stoch(0)=.true.
 122         do i=1,maxflag_stoch
 123           flag_stoch(i)=.false.
 124         enddo
 125 #else
 126         write (iout,*)
 127      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 128 #ifdef MPI
 129         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 130 #endif
 131         stop
 132 #endif
 133       else if (lang.eq.1 .or. lang.eq.4) then
 134         call setup_fricmat
 135       endif
 136 #ifdef MPI
 137       t_langsetup=MPI_Wtime()-tt0
 138       tt0=MPI_Wtime()
 139 #else
 140       t_langsetup=tcpu()-tt0
 141       tt0=tcpu()
 142 #endif
 143       do itime=1,n_timestep
 144         if (ovrtim()) exit
 145         if (large.and. mod(itime,ntwe).eq.0)
 146      &    write (iout,*) "itime",itime
 147         rstcount=rstcount+1
 148         if (lang.gt.0 .and. surfarea .and.
 149      &      mod(itime,reset_fricmat).eq.0) then
 150           if (lang.eq.2 .or. lang.eq.3) then
 151 #ifndef LANG0
 152             call setup_fricmat
 153             if (lang.eq.2) then
 154               call sd_verlet_p_setup
 155             else
 156               call sd_verlet_ciccotti_setup
 157             endif
 158             do i=1,dimen
 159               do j=1,dimen
 160                 pfric0_mat(i,j,0)=pfric_mat(i,j)
 161                 afric0_mat(i,j,0)=afric_mat(i,j)
 162                 vfric0_mat(i,j,0)=vfric_mat(i,j)
 163                 prand0_mat(i,j,0)=prand_mat(i,j)
 164                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
 165                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
 166               enddo
 167             enddo
 168             flag_stoch(0)=.true.
 169             do i=1,maxflag_stoch
 170               flag_stoch(i)=.false.
 171             enddo
 172 #endif
 173           else if (lang.eq.1 .or. lang.eq.4) then
 174             call setup_fricmat
 175           endif
 176           write (iout,'(a,i10)')
 177      &      "Friction matrix reset based on surface area, itime",itime
 178         endif
 179         if (reset_vel .and. tbf .and. lang.eq.0
 180      &      .and. mod(itime,count_reset_vel).eq.0) then
 181           call random_vel
 182           write(iout,'(a,f20.2)')
 183      &     "Velocities reset to random values, time",totT
 184           do i=0,2*nres
 185             do j=1,3
 186               d_t_old(j,i)=d_t(j,i)
 187             enddo
 188           enddo
 189         endif
 190         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
 191           call inertia_tensor
 192           call vcm_vel(vcm)
 193           do j=1,3
 194              d_t(j,0)=d_t(j,0)-vcm(j)
 195           enddo
 196           call kinetic(EK)
 197           kinetic_T=2.0d0/(dimen3*Rb)*EK
 198           scalfac=dsqrt(T_bath/kinetic_T)
 199           write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
 200           do i=0,2*nres
 201             do j=1,3
 202               d_t_old(j,i)=scalfac*d_t(j,i)
 203             enddo
 204           enddo
 205         endif
 206         if (lang.ne.4) then
 207           if (RESPA) then
 208 c Time-reversible RESPA algorithm
 209 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
 210             call RESPA_step(itime)
 211           else
 212 c Variable time step algorithm.
 213             call velverlet_step(itime)
 214           endif
 215         else
 216 #ifdef BROWN
 217           call brown_step(itime)
 218 #else
 219           print *,"Brown dynamics not here!"
 220 #ifdef MPI
 221           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 222 #endif
 223           stop
 224 #endif
 225         endif
 226         if (ntwe.ne.0) then
 227          if (mod(itime,ntwe).eq.0) then
 228            call statout(itime)
 229 C           call enerprint(potEcomp)
 230 C           print *,itime,'AFM',Eafmforc,etot
 231          endif
 232 #ifdef VOUT
 233         do j=1,3
 234           v_work(j)=d_t(j,0)
 235         enddo
 236         ind=3
 237         do i=nnt,nct-1
 238           do j=1,3
 239             ind=ind+1
 240             v_work(ind)=d_t(j,i)
 241           enddo
 242         enddo
 243         do i=nnt,nct
 244           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 245             do j=1,3
 246               ind=ind+1
 247               v_work(ind)=d_t(j,i+nres)
 248             enddo
 249           endif
 250         enddo
 251
 252         write (66,'(80f10.5)')
 253      &    ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
 254         do i=1,ind
 255           v_transf(i)=0.0d0
 256           do j=1,ind
 257             v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
 258           enddo
 259            v_transf(i)= v_transf(i)*dsqrt(geigen(i))
 260         enddo
 261         write (67,'(80f10.5)') (v_transf(i),i=1,ind)
 262 #endif
 263         endif
 264         if (mod(itime,ntwx).eq.0) then
 265           write(iout,*) 'time=',itime
 266 C          call check_ecartint
 267           call returnbox
 268           write (tytul,'("time",f8.2)') totT
 269           if(mdpdb) then
 270              call hairpin(.true.,nharp,iharp)
 271              call secondary2(.true.)
 272              call pdbout(potE,tytul,ipdb)
 273           else
 274              call cartout(totT)
 275           endif
 276         endif
 277         if (rstcount.eq.1000.or.itime.eq.n_timestep) then
 278            open(irest2,file=rest2name,status='unknown')
 279            write(irest2,*) totT,EK,potE,totE,t_bath
 280            do i=1,2*nres
 281             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
 282            enddo
 283            do i=1,2*nres
 284             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
 285            enddo
 286           close(irest2)
 287           rstcount=0
 288         endif
 289       enddo
 290 #ifdef MPI
 291       t_MD=MPI_Wtime()-tt0
 292 #else
 293       t_MD=tcpu()-tt0
 294 #endif
 295       write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
 296      &  '  Timing  ',
 297      & 'MD calculations setup:',t_MDsetup,
 298      & 'Energy & gradient evaluation:',t_enegrad,
 299      & 'Stochastic MD setup:',t_langsetup,
 300      & 'Stochastic MD step setup:',t_sdsetup,
 301      & 'MD steps:',t_MD
 302       write (iout,'(/28(1h=),a25,27(1h=))')
 303      & '  End of MD calculation  '
 304 #ifdef TIMING_ENE
 305       write (iout,*) "time for etotal",t_etotal," elong",t_elong,
 306      &  " eshort",t_eshort
 307       write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
 308      & " time_fricmatmult",time_fricmatmult," time_fsample ",
 309      & time_fsample
 310 #endif
 311       return
 312       end
 313 c-------------------------------------------------------------------------------
 314       subroutine velverlet_step(itime)
 315 c-------------------------------------------------------------------------------
 316 c  Perform a single velocity Verlet step; the time step can be rescaled if
 317 c  increments in accelerations exceed the threshold
 318 c-------------------------------------------------------------------------------
 319       implicit real*8 (a-h,o-z)
 320       include 'DIMENSIONS'
 321 #ifdef MPI
 322       include 'mpif.h'
 323       integer ierror,ierrcode
 324 #endif
 325       include 'COMMON.SETUP'
 326       include 'COMMON.CONTROL'
 327       include 'COMMON.VAR'
 328       include 'COMMON.MD'
 329 #ifndef LANG0
 330       include 'COMMON.LANGEVIN'
 331 #else
 332       include 'COMMON.LANGEVIN.lang0'
 333 #endif
 334       include 'COMMON.CHAIN'
 335       include 'COMMON.DERIV'
 336       include 'COMMON.GEO'
 337       include 'COMMON.LOCAL'
 338       include 'COMMON.INTERACT'
 339       include 'COMMON.IOUNITS'
 340       include 'COMMON.NAMES'
 341       include 'COMMON.TIME1'
 342       include 'COMMON.MUCA'
 343       double precision vcm(3),incr(3)
 344       double precision cm(3),L(3)
 345       integer ilen,count,rstcount
 346       external ilen
 347       character*50 tytul
 348       integer maxcount_scale /20/
 349       common /gucio/ cm
 350       double precision stochforcvec(MAXRES6)
 351       common /stochcalc/ stochforcvec
 352       integer itime
 353       logical scale
 354 c
 355       scale=.true.
 356       icount_scale=0
 357       if (lang.eq.1) then
 358         call sddir_precalc
 359       else if (lang.eq.2 .or. lang.eq.3) then
 360 #ifndef LANG0
 361         call stochastic_force(stochforcvec)
 362 #else
 363         write (iout,*)
 364      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 365 #ifdef MPI
 366         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 367 #endif
 368         stop
 369 #endif
 370       endif
 371       itime_scal=0
 372       do while (scale)
 373         icount_scale=icount_scale+1
 374         if (icount_scale.gt.maxcount_scale) then
 375           write (iout,*)
 376      &      "ERROR: too many attempts at scaling down the time step. ",
 377      &      "amax=",amax,"epdrift=",epdrift,
 378      &      "damax=",damax,"edriftmax=",edriftmax,
 379      &      "d_time=",d_time
 380           call flush(iout)
 381 #ifdef MPI
 382           call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
 383 #endif
 384           stop
 385         endif
 386 c First step of the velocity Verlet algorithm
 387         if (lang.eq.2) then
 388 #ifndef LANG0
 389           call sd_verlet1
 390 #endif
 391         else if (lang.eq.3) then
 392 #ifndef LANG0
 393           call sd_verlet1_ciccotti
 394 #endif
 395         else if (lang.eq.1) then
 396           call sddir_verlet1
 397         else
 398           call verlet1
 399         endif
 400 c Build the chain from the newly calculated coordinates
 401         call chainbuild_cart
 402         if (rattle) call rattle1
 403         if (ntwe.ne.0) then
 404         if (large.and. mod(itime,ntwe).eq.0) then
 405           write (iout,*) "Cartesian and internal coordinates: step 1"
 406           call cartprint
 407           call intout
 408           write (iout,*) "dC"
 409           do i=0,nres
 410             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
 411      &      (dc(j,i+nres),j=1,3)
 412           enddo
 413           write (iout,*) "Accelerations"
 414           do i=0,nres
 415             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 416      &      (d_a(j,i+nres),j=1,3)
 417           enddo
 418           write (iout,*) "Velocities, step 1"
 419           do i=0,nres
 420             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 421      &      (d_t(j,i+nres),j=1,3)
 422           enddo
 423         endif
 424         endif
 425 #ifdef MPI
 426         tt0 = MPI_Wtime()
 427 #else
 428         tt0 = tcpu()
 429 #endif
 430 c Calculate energy and forces
 431         call zerograd
 432         call etotal(potEcomp)
 433         if (large.and. mod(itime,ntwe).eq.0)
 434      &    call enerprint(potEcomp)
 435 #ifdef TIMING_ENE
 436 #ifdef MPI
 437         t_etotal=t_etotal+MPI_Wtime()-tt0
 438 #else
 439         t_etotal=t_etotal+tcpu()-tt0
 440 #endif
 441 #endif
 442         potE=potEcomp(0)-potEcomp(20)
 443         call cartgrad
 444 c Get the new accelerations
 445         call lagrangian
 446 #ifdef MPI
 447         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 448 #else
 449         t_enegrad=t_enegrad+tcpu()-tt0
 450 #endif
 451 c Determine maximum acceleration and scale down the timestep if needed
 452         call max_accel
 453         amax=amax/(itime_scal+1)**2
 454         call predict_edrift(epdrift)
 455         if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
 456 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
 457           scale=.true.
 458           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
 459      &      /dlog(2.0d0)+1
 460           itime_scal=itime_scal+ifac_time
 461 c          fac_time=dmin1(damax/amax,0.5d0)
 462           fac_time=0.5d0**ifac_time
 463           d_time=d_time*fac_time
 464           if (lang.eq.2 .or. lang.eq.3) then
 465 #ifndef LANG0
 466 c            write (iout,*) "Calling sd_verlet_setup: 1"
 467 c Rescale the stochastic forces and recalculate or restore
 468 c the matrices of tinker integrator
 469             if (itime_scal.gt.maxflag_stoch) then
 470               if (large) write (iout,'(a,i5,a)')
 471      &         "Calculate matrices for stochastic step;",
 472      &         " itime_scal ",itime_scal
 473               if (lang.eq.2) then
 474                 call sd_verlet_p_setup
 475               else
 476                 call sd_verlet_ciccotti_setup
 477               endif
 478               write (iout,'(2a,i3,a,i3,1h.)')
 479      &         "Warning: cannot store matrices for stochastic",
 480      &         " integration because the index",itime_scal,
 481      &         " is greater than",maxflag_stoch
 482               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
 483      &         " integration Langevin algorithm for better efficiency."
 484             else if (flag_stoch(itime_scal)) then
 485               if (large) write (iout,'(a,i5,a,l1)')
 486      &         "Restore matrices for stochastic step; itime_scal ",
 487      &         itime_scal," flag ",flag_stoch(itime_scal)
 488               do i=1,dimen
 489                 do j=1,dimen
 490                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
 491                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
 492                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
 493                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
 494                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
 495                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
 496                 enddo
 497               enddo
 498             else
 499               if (large) write (iout,'(2a,i5,a,l1)')
 500      &         "Calculate & store matrices for stochastic step;",
 501      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
 502               if (lang.eq.2) then
 503                 call sd_verlet_p_setup
 504               else
 505                 call sd_verlet_ciccotti_setup
 506               endif
 507               flag_stoch(ifac_time)=.true.
 508               do i=1,dimen
 509                 do j=1,dimen
 510                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
 511                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
 512                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
 513                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
 514                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
 515                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
 516                 enddo
 517               enddo
 518             endif
 519             fac_time=1.0d0/dsqrt(fac_time)
 520             do i=1,dimen
 521               stochforcvec(i)=fac_time*stochforcvec(i)
 522             enddo
 523 #endif
 524           else if (lang.eq.1) then
 525 c Rescale the accelerations due to stochastic forces
 526             fac_time=1.0d0/dsqrt(fac_time)
 527             do i=1,dimen
 528               d_as_work(i)=d_as_work(i)*fac_time
 529             enddo
 530           endif
 531           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
 532      &      "itime",itime," Timestep scaled down to ",
 533      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
 534         else
 535 c Second step of the velocity Verlet algorithm
 536           if (lang.eq.2) then
 537 #ifndef LANG0
 538             call sd_verlet2
 539 #endif
 540           else if (lang.eq.3) then
 541 #ifndef LANG0
 542             call sd_verlet2_ciccotti
 543 #endif
 544           else if (lang.eq.1) then
 545             call sddir_verlet2
 546           else
 547             call verlet2
 548           endif
 549           if (rattle) call rattle2
 550           totT=totT+d_time
 551           totTafm=totT
 552 C          print *,totTafm,"TU?"
 553           if (d_time.ne.d_time0) then
 554             d_time=d_time0
 555 #ifndef   LANG0
 556             if (lang.eq.2 .or. lang.eq.3) then
 557               if (large) write (iout,'(a)')
 558      &         "Restore original matrices for stochastic step"
 559 c              write (iout,*) "Calling sd_verlet_setup: 2"
 560 c Restore the matrices of tinker integrator if the time step has been restored
 561               do i=1,dimen
 562                 do j=1,dimen
 563                   pfric_mat(i,j)=pfric0_mat(i,j,0)
 564                   afric_mat(i,j)=afric0_mat(i,j,0)
 565                   vfric_mat(i,j)=vfric0_mat(i,j,0)
 566                   prand_mat(i,j)=prand0_mat(i,j,0)
 567                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
 568                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
 569                 enddo
 570               enddo
 571             endif
 572 #endif
 573           endif
 574           scale=.false.
 575         endif
 576       enddo
 577 c Calculate the kinetic and the total energy and the kinetic temperature
 578       call kinetic(EK)
 579       totE=EK+potE
 580 c diagnostics
 581 c      call kinetic1(EK1)
 582 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
 583 c end diagnostics
 584 c Couple the system to Berendsen bath if needed
 585       if (tbf .and. lang.eq.0) then
 586         call verlet_bath
 587       endif
 588       kinetic_T=2.0d0/(dimen3*Rb)*EK
 589 c Backup the coordinates, velocities, and accelerations
 590       do i=0,2*nres
 591         do j=1,3
 592           dc_old(j,i)=dc(j,i)
 593           d_t_old(j,i)=d_t(j,i)
 594           d_a_old(j,i)=d_a(j,i)
 595         enddo
 596       enddo
 597       if (ntwe.ne.0) then
 598       if (mod(itime,ntwe).eq.0 .and. large) then
 599         write (iout,*) "Velocities, step 2"
 600         do i=0,nres
 601           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 602      &      (d_t(j,i+nres),j=1,3)
 603         enddo
 604       endif
 605       endif
 606       return
 607       end
 608 c-------------------------------------------------------------------------------
 609       subroutine RESPA_step(itime)
 610 c-------------------------------------------------------------------------------
 611 c  Perform a single RESPA step.
 612 c-------------------------------------------------------------------------------
 613       implicit real*8 (a-h,o-z)
 614       include 'DIMENSIONS'
 615 #ifdef MPI
 616       include 'mpif.h'
 617       integer IERROR,ERRCODE
 618 #endif
 619       include 'COMMON.SETUP'
 620       include 'COMMON.CONTROL'
 621       include 'COMMON.VAR'
 622       include 'COMMON.MD'
 623 #ifndef LANG0
 624       include 'COMMON.LANGEVIN'
 625 #else
 626       include 'COMMON.LANGEVIN.lang0'
 627 #endif
 628       include 'COMMON.CHAIN'
 629       include 'COMMON.DERIV'
 630       include 'COMMON.GEO'
 631       include 'COMMON.LOCAL'
 632       include 'COMMON.INTERACT'
 633       include 'COMMON.IOUNITS'
 634       include 'COMMON.NAMES'
 635       include 'COMMON.TIME1'
 636       double precision energia_short(0:n_ene),
 637      & energia_long(0:n_ene)
 638       double precision cm(3),L(3),vcm(3),incr(3)
 639       double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
 640      & d_a_old0(3,0:maxres2)
 641       logical PRINT_AMTS_MSG /.false./
 642       integer ilen,count,rstcount
 643       external ilen
 644       character*50 tytul
 645       integer maxcount_scale /10/
 646       common /gucio/ cm
 647       double precision stochforcvec(MAXRES6)
 648       common /stochcalc/ stochforcvec
 649       integer itime
 650       logical scale
 651       common /cipiszcze/ itt
 652       itt=itime
 653       if (ntwe.ne.0) then
 654       if (large.and. mod(itime,ntwe).eq.0) then
 655         write (iout,*) "***************** RESPA itime",itime
 656         write (iout,*) "Cartesian and internal coordinates: step 0"
 657 c        call cartprint
 658         call pdbout(0.0d0,"cipiszcze",iout)
 659         call intout
 660         write (iout,*) "Accelerations from long-range forces"
 661         do i=0,nres
 662           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 663      &      (d_a(j,i+nres),j=1,3)
 664         enddo
 665         write (iout,*) "Velocities, step 0"
 666         do i=0,nres
 667           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 668      &      (d_t(j,i+nres),j=1,3)
 669         enddo
 670       endif
 671       endif
 672 c
 673 c Perform the initial RESPA step (increment velocities)
 674 c      write (iout,*) "*********************** RESPA ini"
 675       call RESPA_vel
 676       if (ntwe.ne.0) then
 677       if (mod(itime,ntwe).eq.0 .and. large) then
 678         write (iout,*) "Velocities, end"
 679         do i=0,nres
 680           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 681      &      (d_t(j,i+nres),j=1,3)
 682         enddo
 683       endif
 684       endif
 685 c Compute the short-range forces
 686 #ifdef MPI
 687       tt0 =MPI_Wtime()
 688 #else
 689       tt0 = tcpu()
 690 #endif
 691 C 7/2/2009 commented out
 692 c      call zerograd
 693 c      call etotal_short(energia_short)
 694 c      call cartgrad
 695 c      call lagrangian
 696 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
 697         do i=0,2*nres
 698           do j=1,3
 699             d_a(j,i)=d_a_short(j,i)
 700           enddo
 701         enddo
 702       if (ntwe.ne.0) then
 703       if (large.and. mod(itime,ntwe).eq.0) then
 704         write (iout,*) "energia_short",energia_short(0)
 705         write (iout,*) "Accelerations from short-range forces"
 706         do i=0,nres
 707           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 708      &      (d_a(j,i+nres),j=1,3)
 709         enddo
 710       endif
 711       endif
 712 #ifdef MPI
 713         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 714 #else
 715         t_enegrad=t_enegrad+tcpu()-tt0
 716 #endif
 717       do i=0,2*nres
 718         do j=1,3
 719           dc_old(j,i)=dc(j,i)
 720           d_t_old(j,i)=d_t(j,i)
 721           d_a_old(j,i)=d_a(j,i)
 722         enddo
 723       enddo
 724 c 6/30/08 A-MTS: attempt at increasing the split number
 725       do i=0,2*nres
 726         do j=1,3
 727           dc_old0(j,i)=dc_old(j,i)
 728           d_t_old0(j,i)=d_t_old(j,i)
 729           d_a_old0(j,i)=d_a_old(j,i)
 730         enddo
 731       enddo
 732       if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
 733       if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
 734 c
 735       scale=.true.
 736       d_time0=d_time
 737       do while (scale)
 738
 739       scale=.false.
 740 c      write (iout,*) "itime",itime," ntime_split",ntime_split
 741 c Split the time step
 742       d_time=d_time0/ntime_split
 743 c Perform the short-range RESPA steps (velocity Verlet increments of
 744 c positions and velocities using short-range forces)
 745 c      write (iout,*) "*********************** RESPA split"
 746       do itsplit=1,ntime_split
 747         if (lang.eq.1) then
 748           call sddir_precalc
 749         else if (lang.eq.2 .or. lang.eq.3) then
 750 #ifndef LANG0
 751           call stochastic_force(stochforcvec)
 752 #else
 753           write (iout,*)
 754      &      "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 755 #ifdef MPI
 756           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 757 #endif
 758           stop
 759 #endif
 760         endif
 761 c First step of the velocity Verlet algorithm
 762         if (lang.eq.2) then
 763 #ifndef LANG0
 764           call sd_verlet1
 765 #endif
 766         else if (lang.eq.3) then
 767 #ifndef LANG0
 768           call sd_verlet1_ciccotti
 769 #endif
 770         else if (lang.eq.1) then
 771           call sddir_verlet1
 772         else
 773           call verlet1
 774         endif
 775 c Build the chain from the newly calculated coordinates
 776         call chainbuild_cart
 777         if (rattle) call rattle1
 778         if (ntwe.ne.0) then
 779         if (large.and. mod(itime,ntwe).eq.0) then
 780           write (iout,*) "***** ITSPLIT",itsplit
 781           write (iout,*) "Cartesian and internal coordinates: step 1"
 782           call pdbout(0.0d0,"cipiszcze",iout)
 783 c          call cartprint
 784           call intout
 785           write (iout,*) "Velocities, step 1"
 786           do i=0,nres
 787             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 788      &        (d_t(j,i+nres),j=1,3)
 789           enddo
 790         endif
 791         endif
 792 #ifdef MPI
 793         tt0 = MPI_Wtime()
 794 #else
 795         tt0 = tcpu()
 796 #endif
 797 c Calculate energy and forces
 798         call zerograd
 799         call etotal_short(energia_short)
 800         if (large.and. mod(itime,ntwe).eq.0)
 801      &    call enerprint(energia_short)
 802 #ifdef TIMING_ENE
 803 #ifdef MPI
 804         t_eshort=t_eshort+MPI_Wtime()-tt0
 805 #else
 806         t_eshort=t_eshort+tcpu()-tt0
 807 #endif
 808 #endif
 809         call cartgrad
 810 c Get the new accelerations
 811         call lagrangian
 812 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
 813         do i=0,2*nres
 814           do j=1,3
 815             d_a_short(j,i)=d_a(j,i)
 816           enddo
 817         enddo
 818         if (ntwe.ne.0) then
 819         if (large.and. mod(itime,ntwe).eq.0) then
 820           write (iout,*)"energia_short",energia_short(0)
 821           write (iout,*) "Accelerations from short-range forces"
 822           do i=0,nres
 823             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 824      &        (d_a(j,i+nres),j=1,3)
 825           enddo
 826         endif
 827         endif
 828 c 6/30/08 A-MTS
 829 c Determine maximum acceleration and scale down the timestep if needed
 830         call max_accel
 831         amax=amax/ntime_split**2
 832         call predict_edrift(epdrift)
 833         if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
 834      &   write (iout,*) "amax",amax," damax",damax,
 835      &   " epdrift",epdrift," epdriftmax",epdriftmax
 836 c Exit loop and try with increased split number if the change of
 837 c acceleration is too big
 838         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
 839           if (ntime_split.lt.maxtime_split) then
 840             scale=.true.
 841             ntime_split=ntime_split*2
 842             do i=0,2*nres
 843               do j=1,3
 844                 dc_old(j,i)=dc_old0(j,i)
 845                 d_t_old(j,i)=d_t_old0(j,i)
 846                 d_a_old(j,i)=d_a_old0(j,i)
 847               enddo
 848             enddo
 849             if (PRINT_AMTS_MSG) then
 850             write (iout,*) "acceleration/energy drift too large",amax,
 851      &      epdrift," split increased to ",ntime_split," itime",itime,
 852      &       " itsplit",itsplit
 853             endif
 854             exit
 855           else
 856             write (iout,*)
 857      &      "Uh-hu. Bumpy landscape. Maximum splitting number",
 858      &       maxtime_split,
 859      &      " already reached!!! Trying to carry on!"
 860           endif
 861         endif
 862 #ifdef MPI
 863         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 864 #else
 865         t_enegrad=t_enegrad+tcpu()-tt0
 866 #endif
 867 c Second step of the velocity Verlet algorithm
 868         if (lang.eq.2) then
 869 #ifndef LANG0
 870           call sd_verlet2
 871 #endif
 872         else if (lang.eq.3) then
 873 #ifndef LANG0
 874           call sd_verlet2_ciccotti
 875 #endif
 876         else if (lang.eq.1) then
 877           call sddir_verlet2
 878         else
 879           call verlet2
 880         endif
 881         if (rattle) call rattle2
 882 c Backup the coordinates, velocities, and accelerations
 883         do i=0,2*nres
 884           do j=1,3
 885             dc_old(j,i)=dc(j,i)
 886             d_t_old(j,i)=d_t(j,i)
 887             d_a_old(j,i)=d_a(j,i)
 888           enddo
 889         enddo
 890       enddo
 891
 892       enddo ! while scale
 893
 894 c Restore the time step
 895       d_time=d_time0
 896 c Compute long-range forces
 897 #ifdef MPI
 898       tt0 =MPI_Wtime()
 899 #else
 900       tt0 = tcpu()
 901 #endif
 902       call zerograd
 903       call etotal_long(energia_long)
 904       if (large.and. mod(itime,ntwe).eq.0)
 905      &    call enerprint(energia_long)
 906 #ifdef TIMING_ENE
 907 #ifdef MPI
 908         t_elong=t_elong+MPI_Wtime()-tt0
 909 #else
 910         t_elong=t_elong+tcpu()-tt0
 911 #endif
 912 #endif
 913       call cartgrad
 914       call lagrangian
 915 #ifdef MPI
 916         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 917 #else
 918         t_enegrad=t_enegrad+tcpu()-tt0
 919 #endif
 920 c Compute accelerations from long-range forces
 921       if (ntwe.ne.0) then
 922       if (large.and. mod(itime,ntwe).eq.0) then
 923         write (iout,*) "energia_long",energia_long(0)
 924         write (iout,*) "Cartesian and internal coordinates: step 2"
 925 c        call cartprint
 926         call pdbout(0.0d0,"cipiszcze",iout)
 927         call intout
 928         write (iout,*) "Accelerations from long-range forces"
 929         do i=0,nres
 930           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 931      &      (d_a(j,i+nres),j=1,3)
 932         enddo
 933         write (iout,*) "Velocities, step 2"
 934         do i=0,nres
 935           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 936      &      (d_t(j,i+nres),j=1,3)
 937         enddo
 938       endif
 939       endif
 940 c Compute the final RESPA step (increment velocities)
 941 c      write (iout,*) "*********************** RESPA fin"
 942       call RESPA_vel
 943 c Compute the complete potential energy
 944       do i=0,n_ene
 945         potEcomp(i)=energia_short(i)+energia_long(i)
 946       enddo
 947       potE=potEcomp(0)-potEcomp(20)
 948 c      potE=energia_short(0)+energia_long(0)
 949       totT=totT+d_time
 950       totTafm=totT
 951 c Calculate the kinetic and the total energy and the kinetic temperature
 952       call kinetic(EK)
 953       totE=EK+potE
 954 c Couple the system to Berendsen bath if needed
 955       if (tbf .and. lang.eq.0) then
 956         call verlet_bath
 957       endif
 958       kinetic_T=2.0d0/(dimen3*Rb)*EK
 959 c Backup the coordinates, velocities, and accelerations
 960       if (ntwe.ne.0) then
 961       if (mod(itime,ntwe).eq.0 .and. large) then
 962         write (iout,*) "Velocities, end"
 963         do i=0,nres
 964           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 965      &      (d_t(j,i+nres),j=1,3)
 966         enddo
 967       endif
 968       endif
 969       return
 970       end
 971 c---------------------------------------------------------------------
 972       subroutine RESPA_vel
 973 c  First and last RESPA step (incrementing velocities using long-range
 974 c  forces).
 975       implicit real*8 (a-h,o-z)
 976       include 'DIMENSIONS'
 977       include 'COMMON.CONTROL'
 978       include 'COMMON.VAR'
 979       include 'COMMON.MD'
 980       include 'COMMON.CHAIN'
 981       include 'COMMON.DERIV'
 982       include 'COMMON.GEO'
 983       include 'COMMON.LOCAL'
 984       include 'COMMON.INTERACT'
 985       include 'COMMON.IOUNITS'
 986       include 'COMMON.NAMES'
 987       do j=1,3
 988         d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
 989       enddo
 990       do i=nnt,nct-1
 991         do j=1,3
 992           d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
 993         enddo
 994       enddo
 995       do i=nnt,nct
 996         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 997           inres=i+nres
 998           do j=1,3
 999             d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1000           enddo
1001         endif
1002       enddo
1003       return
1004       end
1005 c-----------------------------------------------------------------
1006       subroutine verlet1
1007 c Applying velocity Verlet algorithm - step 1 to coordinates
1008       implicit real*8 (a-h,o-z)
1009       include 'DIMENSIONS'
1010       include 'COMMON.CONTROL'
1011       include 'COMMON.VAR'
1012       include 'COMMON.MD'
1013       include 'COMMON.CHAIN'
1014       include 'COMMON.DERIV'
1015       include 'COMMON.GEO'
1016       include 'COMMON.LOCAL'
1017       include 'COMMON.INTERACT'
1018       include 'COMMON.IOUNITS'
1019       include 'COMMON.NAMES'
1020       double precision adt,adt2
1021
1022 #ifdef DEBUG
1023       write (iout,*) "VELVERLET1 START: DC"
1024       do i=0,nres
1025         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1026      &   (dc(j,i+nres),j=1,3)
1027       enddo
1028 #endif
1029       do j=1,3
1030         adt=d_a_old(j,0)*d_time
1031         adt2=0.5d0*adt
1032         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1033         d_t_new(j,0)=d_t_old(j,0)+adt2
1034         d_t(j,0)=d_t_old(j,0)+adt
1035       enddo
1036       do i=nnt,nct-1
1037 C      SPYTAC ADAMA
1038 C       do i=0,nres
1039         do j=1,3
1040           adt=d_a_old(j,i)*d_time
1041           adt2=0.5d0*adt
1042           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1043           d_t_new(j,i)=d_t_old(j,i)+adt2
1044           d_t(j,i)=d_t_old(j,i)+adt
1045         enddo
1046       enddo
1047       do i=nnt,nct
1048 C        do i=0,nres
1049         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1050           inres=i+nres
1051           do j=1,3
1052             adt=d_a_old(j,inres)*d_time
1053             adt2=0.5d0*adt
1054             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1055             d_t_new(j,inres)=d_t_old(j,inres)+adt2
1056             d_t(j,inres)=d_t_old(j,inres)+adt
1057           enddo
1058         endif
1059       enddo
1060 #ifdef DEBUG
1061       write (iout,*) "VELVERLET1 END: DC"
1062       do i=0,nres
1063         write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1064      &   (dc(j,i+nres),j=1,3)
1065       enddo
1066 #endif
1067       return
1068       end
1069 c---------------------------------------------------------------------
1070       subroutine verlet2
1071 c  Step 2 of the velocity Verlet algorithm: update velocities
1072       implicit real*8 (a-h,o-z)
1073       include 'DIMENSIONS'
1074       include 'COMMON.CONTROL'
1075       include 'COMMON.VAR'
1076       include 'COMMON.MD'
1077       include 'COMMON.CHAIN'
1078       include 'COMMON.DERIV'
1079       include 'COMMON.GEO'
1080       include 'COMMON.LOCAL'
1081       include 'COMMON.INTERACT'
1082       include 'COMMON.IOUNITS'
1083       include 'COMMON.NAMES'
1084       do j=1,3
1085         d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1086       enddo
1087       do i=nnt,nct-1
1088         do j=1,3
1089           d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1090         enddo
1091       enddo
1092       do i=nnt,nct
1093         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1094           inres=i+nres
1095           do j=1,3
1096             d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1097           enddo
1098         endif
1099       enddo
1100       return
1101       end
1102 c-----------------------------------------------------------------
1103       subroutine sddir_precalc
1104 c Applying velocity Verlet algorithm - step 1 to coordinates
1105       implicit real*8 (a-h,o-z)
1106       include 'DIMENSIONS'
1107 #ifdef MPI
1108       include 'mpif.h'
1109 #endif
1110       include 'COMMON.CONTROL'
1111       include 'COMMON.VAR'
1112       include 'COMMON.MD'
1113 #ifndef LANG0
1114       include 'COMMON.LANGEVIN'
1115 #else
1116       include 'COMMON.LANGEVIN.lang0'
1117 #endif
1118       include 'COMMON.CHAIN'
1119       include 'COMMON.DERIV'
1120       include 'COMMON.GEO'
1121       include 'COMMON.LOCAL'
1122       include 'COMMON.INTERACT'
1123       include 'COMMON.IOUNITS'
1124       include 'COMMON.NAMES'
1125       include 'COMMON.TIME1'
1126       double precision stochforcvec(MAXRES6)
1127       common /stochcalc/ stochforcvec
1128 c
1129 c Compute friction and stochastic forces
1130 c
1131 #ifdef MPI
1132       time00=MPI_Wtime()
1133 #else
1134       time00=tcpu()
1135 #endif
1136       call friction_force
1137 #ifdef MPI
1138       time_fric=time_fric+MPI_Wtime()-time00
1139       time00=MPI_Wtime()
1140 #else
1141       time_fric=time_fric+tcpu()-time00
1142       time00=tcpu()
1143 #endif
1144       call stochastic_force(stochforcvec)
1145 #ifdef MPI
1146       time_stoch=time_stoch+MPI_Wtime()-time00
1147 #else
1148       time_stoch=time_stoch+tcpu()-time00
1149 #endif
1150 c
1151 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1152 c forces (d_as_work)
1153 c
1154       call ginv_mult(fric_work, d_af_work)
1155       call ginv_mult(stochforcvec, d_as_work)
1156       return
1157       end
1158 c---------------------------------------------------------------------
1159       subroutine sddir_verlet1
1160 c Applying velocity Verlet algorithm - step 1 to velocities
1161       implicit real*8 (a-h,o-z)
1162       include 'DIMENSIONS'
1163       include 'COMMON.CONTROL'
1164       include 'COMMON.VAR'
1165       include 'COMMON.MD'
1166 #ifndef LANG0
1167       include 'COMMON.LANGEVIN'
1168 #else
1169       include 'COMMON.LANGEVIN.lang0'
1170 #endif
1171       include 'COMMON.CHAIN'
1172       include 'COMMON.DERIV'
1173       include 'COMMON.GEO'
1174       include 'COMMON.LOCAL'
1175       include 'COMMON.INTERACT'
1176       include 'COMMON.IOUNITS'
1177       include 'COMMON.NAMES'
1178 c Revised 3/31/05 AL: correlation between random contributions to
1179 c position and velocity increments included.
1180       double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
1181       double precision adt,adt2
1182 c
1183 c Add the contribution from BOTH friction and stochastic force to the
1184 c coordinates, but ONLY the contribution from the friction forces to velocities
1185 c
1186       do j=1,3
1187         adt=(d_a_old(j,0)+d_af_work(j))*d_time
1188         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1189         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1190         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1191         d_t(j,0)=d_t_old(j,0)+adt
1192       enddo
1193       ind=3
1194       do i=nnt,nct-1
1195         do j=1,3
1196           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1197           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1198           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1199           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1200           d_t(j,i)=d_t_old(j,i)+adt
1201         enddo
1202         ind=ind+3
1203       enddo
1204       do i=nnt,nct
1205         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1206           inres=i+nres
1207           do j=1,3
1208             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1209             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1210             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1211             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1212             d_t(j,inres)=d_t_old(j,inres)+adt
1213           enddo
1214           ind=ind+3
1215         endif
1216       enddo
1217       return
1218       end
1219 c---------------------------------------------------------------------
1220       subroutine sddir_verlet2
1221 c  Calculating the adjusted velocities for accelerations
1222       implicit real*8 (a-h,o-z)
1223       include 'DIMENSIONS'
1224       include 'COMMON.CONTROL'
1225       include 'COMMON.VAR'
1226       include 'COMMON.MD'
1227 #ifndef LANG0
1228       include 'COMMON.LANGEVIN'
1229 #else
1230       include 'COMMON.LANGEVIN.lang0'
1231 #endif
1232       include 'COMMON.CHAIN'
1233       include 'COMMON.DERIV'
1234       include 'COMMON.GEO'
1235       include 'COMMON.LOCAL'
1236       include 'COMMON.INTERACT'
1237       include 'COMMON.IOUNITS'
1238       include 'COMMON.NAMES'
1239       double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
1240       double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
1241 c Revised 3/31/05 AL: correlation between random contributions to
1242 c position and velocity increments included.
1243 c The correlation coefficients are calculated at low-friction limit.
1244 c Also, friction forces are now not calculated with new velocities.
1245
1246 c      call friction_force
1247       call stochastic_force(stochforcvec)
1248 c
1249 c Compute the acceleration due to friction forces (d_af_work) and stochastic
1250 c forces (d_as_work)
1251 c
1252       call ginv_mult(stochforcvec, d_as_work1)
1253
1254 c
1255 c Update velocities
1256 c
1257       do j=1,3
1258         d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
1259      &    +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1260       enddo
1261       ind=3
1262       do i=nnt,nct-1
1263         do j=1,3
1264           d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
1265      &     +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
1266         enddo
1267         ind=ind+3
1268       enddo
1269       do i=nnt,nct
1270         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1271           inres=i+nres
1272           do j=1,3
1273             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
1274      &       +d_af_work(ind+j))+sin60*d_as_work(ind+j)
1275      &       +cos60*d_as_work1(ind+j))*d_time
1276           enddo
1277           ind=ind+3
1278         endif
1279       enddo
1280       return
1281       end
1282 c---------------------------------------------------------------------
1283       subroutine max_accel
1284 c
1285 c Find the maximum difference in the accelerations of the the sites
1286 c at the beginning and the end of the time step.
1287 c
1288       implicit real*8 (a-h,o-z)
1289       include 'DIMENSIONS'
1290       include 'COMMON.CONTROL'
1291       include 'COMMON.VAR'
1292       include 'COMMON.MD'
1293       include 'COMMON.CHAIN'
1294       include 'COMMON.DERIV'
1295       include 'COMMON.GEO'
1296       include 'COMMON.LOCAL'
1297       include 'COMMON.INTERACT'
1298       include 'COMMON.IOUNITS'
1299       double precision aux(3),accel(3),accel_old(3),dacc
1300       do j=1,3
1301 c        aux(j)=d_a(j,0)-d_a_old(j,0)
1302          accel_old(j)=d_a_old(j,0)
1303          accel(j)=d_a(j,0)
1304       enddo
1305       amax=0.0d0
1306       do i=nnt,nct
1307 c Backbone
1308         if (i.lt.nct) then
1309 c 7/3/08 changed to asymmetric difference
1310           do j=1,3
1311 c            accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
1312             accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
1313             accel(j)=accel(j)+0.5d0*d_a(j,i)
1314 c            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1315             if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1316               dacc=dabs(accel(j)-accel_old(j))
1317 c              write (iout,*) i,dacc
1318               if (dacc.gt.amax) amax=dacc
1319             endif
1320           enddo
1321         endif
1322       enddo
1323 c Side chains
1324       do j=1,3
1325 c        accel(j)=aux(j)
1326         accel_old(j)=d_a_old(j,0)
1327         accel(j)=d_a(j,0)
1328       enddo
1329       if (nnt.eq.2) then
1330         do j=1,3
1331           accel_old(j)=accel_old(j)+d_a_old(j,1)
1332           accel(j)=accel(j)+d_a(j,1)
1333         enddo
1334       endif
1335       do i=nnt,nct
1336         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1337           do j=1,3
1338 c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
1339             accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
1340             accel(j)=accel(j)+d_a(j,i+nres)
1341           enddo
1342         endif
1343         do j=1,3
1344 c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
1345           if (dabs(accel(j)).gt.dabs(accel_old(j))) then
1346             dacc=dabs(accel(j)-accel_old(j))
1347 c            write (iout,*) "side-chain",i,dacc
1348             if (dacc.gt.amax) amax=dacc
1349           endif
1350         enddo
1351         do j=1,3
1352           accel_old(j)=accel_old(j)+d_a_old(j,i)
1353           accel(j)=accel(j)+d_a(j,i)
1354 c          aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
1355         enddo
1356       enddo
1357       return
1358       end
1359 c---------------------------------------------------------------------
1360       subroutine predict_edrift(epdrift)
1361 c
1362 c Predict the drift of the potential energy
1363 c
1364       implicit real*8 (a-h,o-z)
1365       include 'DIMENSIONS'
1366       include 'COMMON.CONTROL'
1367       include 'COMMON.VAR'
1368       include 'COMMON.MD'
1369       include 'COMMON.CHAIN'
1370       include 'COMMON.DERIV'
1371       include 'COMMON.GEO'
1372       include 'COMMON.LOCAL'
1373       include 'COMMON.INTERACT'
1374       include 'COMMON.IOUNITS'
1375       include 'COMMON.MUCA'
1376       double precision epdrift,epdriftij
1377 c Drift of the potential energy
1378       epdrift=0.0d0
1379       do i=nnt,nct
1380 c Backbone
1381         if (i.lt.nct) then
1382           do j=1,3
1383             epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
1384             if (lmuca) epdriftij=epdriftij*factor
1385 c            write (iout,*) "back",i,j,epdriftij
1386             if (epdriftij.gt.epdrift) epdrift=epdriftij
1387           enddo
1388         endif
1389 c Side chains
1390         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1391           do j=1,3
1392             epdriftij=
1393      &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
1394             if (lmuca) epdriftij=epdriftij*factor
1395 c            write (iout,*) "side",i,j,epdriftij
1396             if (epdriftij.gt.epdrift) epdrift=epdriftij
1397           enddo
1398         endif
1399       enddo
1400       epdrift=0.5d0*epdrift*d_time*d_time
1401 c      write (iout,*) "epdrift",epdrift
1402       return
1403       end
1404 c-----------------------------------------------------------------------
1405       subroutine verlet_bath
1406 c
1407 c  Coupling to the thermostat by using the Berendsen algorithm
1408 c
1409       implicit real*8 (a-h,o-z)
1410       include 'DIMENSIONS'
1411       include 'COMMON.CONTROL'
1412       include 'COMMON.VAR'
1413       include 'COMMON.MD'
1414       include 'COMMON.CHAIN'
1415       include 'COMMON.DERIV'
1416       include 'COMMON.GEO'
1417       include 'COMMON.LOCAL'
1418       include 'COMMON.INTERACT'
1419       include 'COMMON.IOUNITS'
1420       include 'COMMON.NAMES'
1421       double precision T_half,fact
1422 c
1423       T_half=2.0d0/(dimen3*Rb)*EK
1424       fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
1425 c      write(iout,*) "T_half", T_half
1426 c      write(iout,*) "EK", EK
1427 c      write(iout,*) "fact", fact
1428       do j=1,3
1429         d_t(j,0)=fact*d_t(j,0)
1430       enddo
1431       do i=nnt,nct-1
1432         do j=1,3
1433           d_t(j,i)=fact*d_t(j,i)
1434         enddo
1435       enddo
1436       do i=nnt,nct
1437         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1438           inres=i+nres
1439           do j=1,3
1440             d_t(j,inres)=fact*d_t(j,inres)
1441           enddo
1442         endif
1443       enddo
1444       return
1445       end
1446 c---------------------------------------------------------
1447       subroutine init_MD
1448 c  Set up the initial conditions of a MD simulation
1449       implicit real*8 (a-h,o-z)
1450       include 'DIMENSIONS'
1451 #ifdef MP
1452       include 'mpif.h'
1453       character*16 form
1454       integer IERROR,ERRCODE
1455 #endif
1456       include 'COMMON.SETUP'
1457       include 'COMMON.CONTROL'
1458       include 'COMMON.VAR'
1459       include 'COMMON.MD'
1460 #ifndef LANG0
1461       include 'COMMON.LANGEVIN'
1462 #else
1463       include 'COMMON.LANGEVIN.lang0'
1464 #endif
1465       include 'COMMON.CHAIN'
1466       include 'COMMON.DERIV'
1467       include 'COMMON.GEO'
1468       include 'COMMON.LOCAL'
1469       include 'COMMON.INTERACT'
1470       include 'COMMON.IOUNITS'
1471       include 'COMMON.NAMES'
1472       include 'COMMON.REMD'
1473       real*8 energia_long(0:n_ene),
1474      &  energia_short(0:n_ene),energia(0:n_ene),vcm(3),incr(3)
1475       double precision cm(3),L(3),xv,sigv,lowb,highb
1476       double precision varia(maxvar)
1477       character*256 qstr
1478       integer ilen
1479       external ilen
1480       character*50 tytul
1481       logical file_exist
1482       common /gucio/ cm
1483       d_time0=d_time
1484 c      write(iout,*) "d_time", d_time
1485 c Compute the standard deviations of stochastic forces for Langevin dynamics
1486 c if the friction coefficients do not depend on surface area
1487       if (lang.gt.0 .and. .not.surfarea) then
1488         do i=nnt,nct-1
1489           stdforcp(i)=stdfp*dsqrt(gamp)
1490         enddo
1491         do i=nnt,nct
1492           stdforcsc(i)=stdfsc(iabs(itype(i)))
1493      &                *dsqrt(gamsc(iabs(itype(i))))
1494         enddo
1495       endif
1496 c Open the pdb file for snapshotshots
1497 #ifdef MPI
1498       if(mdpdb) then
1499         if (ilen(tmpdir).gt.0)
1500      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1501      &      liczba(:ilen(liczba))//".pdb")
1502         open(ipdb,
1503      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1504      &  //".pdb")
1505       else
1506 #ifdef NOXDR
1507         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1508      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1509      &      liczba(:ilen(liczba))//".x")
1510         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1511      &  //".x"
1512 #else
1513         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1514      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1515      &      liczba(:ilen(liczba))//".cx")
1516         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1517      &  //".cx"
1518 #endif
1519       endif
1520 #else
1521       if(mdpdb) then
1522          if (ilen(tmpdir).gt.0)
1523      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1524          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1525       else
1526          if (ilen(tmpdir).gt.0)
1527      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1528          cartname=prefix(:ilen(prefix))//"_MD.cx"
1529       endif
1530 #endif
1531       if (usampl) then
1532         write (qstr,'(256(1h ))')
1533         ipos=1
1534         do i=1,nfrag
1535           iq = qinfrag(i,iset)*10
1536           iw = wfrag(i,iset)/100
1537           if (iw.gt.0) then
1538             if(me.eq.king.or..not.out1file)
1539      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1540             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1541             ipos=ipos+7
1542           endif
1543         enddo
1544         do i=1,npair
1545           iq = qinpair(i,iset)*10
1546           iw = wpair(i,iset)/100
1547           if (iw.gt.0) then
1548             if(me.eq.king.or..not.out1file)
1549      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1550             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1551             ipos=ipos+7
1552           endif
1553         enddo
1554 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1555 #ifdef NOXDR
1556 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1557 #else
1558 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1559 #endif
1560 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1561       endif
1562       icg=1
1563       if (rest) then
1564        if (restart1file) then
1565          if (me.eq.king)
1566      &     inquire(file=mremd_rst_name,exist=file_exist)
1567 #ifdef MPI
1568            write (*,*) me," Before broadcast: file_exist",file_exist
1569          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1570      &          IERR)
1571          write (*,*) me," After broadcast: file_exist",file_exist
1572 c        inquire(file=mremd_rst_name,exist=file_exist)
1573 #endif
1574         if(me.eq.king.or..not.out1file)
1575      &   write(iout,*) "Initial state read by master and distributed"
1576        else
1577          if (ilen(tmpdir).gt.0)
1578      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1579      &      //liczba(:ilen(liczba))//'.rst')
1580         inquire(file=rest2name,exist=file_exist)
1581        endif
1582        if(file_exist) then
1583          if(.not.restart1file) then
1584            if(me.eq.king.or..not.out1file)
1585      &      write(iout,*) "Initial state will be read from file ",
1586      &      rest2name(:ilen(rest2name))
1587            call readrst
1588          endif
1589          call rescale_weights(t_bath)
1590        else
1591         if(me.eq.king.or..not.out1file)then
1592          if (restart1file) then
1593           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1594      &       " does not exist"
1595          else
1596           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1597      &       " does not exist"
1598          endif
1599          write(iout,*) "Initial velocities randomly generated"
1600         endif
1601         call random_vel
1602         totT=0.0d0
1603         totTafm=totT
1604        endif
1605       else
1606 c Generate initial velocities
1607         if(me.eq.king.or..not.out1file)
1608      &   write(iout,*) "Initial velocities randomly generated"
1609         call random_vel
1610         totT=0.0d0
1611 CtotTafm is the variable for AFM time which eclipsed during
1612         totTafm=totT
1613       endif
1614 c      rest2name = prefix(:ilen(prefix))//'.rst'
1615       if(me.eq.king.or..not.out1file)then
1616        write (iout,*) "Initial velocities"
1617        do i=0,nres
1618          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1619      &   (d_t(j,i+nres),j=1,3)
1620        enddo
1621 c  Zeroing the total angular momentum of the system
1622        write(iout,*) "Calling the zero-angular
1623      & momentum subroutine"
1624       endif
1625       call inertia_tensor
1626 c  Getting the potential energy and forces and velocities and accelerations
1627       call vcm_vel(vcm)
1628 c      write (iout,*) "velocity of the center of the mass:"
1629 c      write (iout,*) (vcm(j),j=1,3)
1630       do j=1,3
1631         d_t(j,0)=d_t(j,0)-vcm(j)
1632       enddo
1633 c Removing the velocity of the center of mass
1634       call vcm_vel(vcm)
1635       if(me.eq.king.or..not.out1file)then
1636        write (iout,*) "vcm right after adjustment:"
1637        write (iout,*) (vcm(j),j=1,3)
1638       endif
1639       if (.not.rest) then
1640          if (start_from_model) then
1641           i_model=iran_num(1,constr_homology)
1642           write (iout,*) 'starting from model ',i_model
1643           do i=1,2*nres
1644            do j=1,3
1645             c(j,i)=chomo(j,i,i_model)
1646            enddo
1647           enddo
1648
1649           call int_from_cart(.true.,.false.)
1650           call sc_loc_geom(.false.)
1651           do i=1,nres-1
1652             do j=1,3
1653              dc(j,i)=c(j,i+1)-c(j,i)
1654              dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1655             enddo
1656           enddo
1657           do i=2,nres-1
1658             do j=1,3
1659              dc(j,i+nres)=c(j,i+nres)-c(j,i)
1660              dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1661             enddo
1662           enddo
1663           call etotal(energia(0))
1664          endif
1665 c         call chainbuild_cart
1666          write (iout,*) "Initial energies"
1667          call enerprint(energia(0))
1668          write (iout,*) "PREMINIM ",preminim
1669          if(iranconf.ne.0 .or. preminim) then
1670           if (overlapsc) then
1671            write (iout,*) 'Calling OVERLAP_SC'
1672            call overlap_sc(fail)
1673           endif
1674
1675           if (searchsc) then
1676            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1677            print *,'SC_move',nft_sc,etot
1678            if(me.eq.king.or..not.out1file)
1679      &      write(iout,*) 'SC_move',nft_sc,etot
1680           endif
1681
1682           if(dccart)then
1683            print *, 'Calling MINIM_DC'
1684            call minim_dc(etot,iretcode,nfun)
1685           else
1686            call geom_to_var(nvar,varia)
1687            print *,'Calling MINIMIZE.'
1688            call minimize(etot,varia,iretcode,nfun)
1689            call var_to_geom(nvar,varia)
1690           endif
1691           if(me.eq.king.or..not.out1file) then
1692              write(iout,*) "Minimized energy is",etot
1693              write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1694              call etotal(potEcomp)
1695              call enerprint(potEcomp)
1696           endif
1697          endif
1698          if(iranconf.ne.0) then
1699 c 8/22/17 AL Loop to produce a low-energy random conformation
1700           DO iranmin=1,10
1701           if (overlapsc) then
1702            if(me.eq.king.or..not.out1file)
1703      &        write (iout,*) 'Calling OVERLAP_SC'
1704            call overlap_sc(fail)
1705           endif
1706
1707           if (searchsc) then
1708            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1709            print *,'SC_move',nft_sc,etot
1710            if(me.eq.king.or..not.out1file)
1711      &      write(iout,*) 'SC_move',nft_sc,etot
1712           endif
1713
1714           if(dccart)then
1715            print *, 'Calling MINIM_DC'
1716            call minim_dc(etot,iretcode,nfun)
1717           else
1718            call geom_to_var(nvar,varia)
1719            print *,'Calling MINIMIZE.'
1720            call minimize(etot,varia,iretcode,nfun)
1721            call var_to_geom(nvar,varia)
1722           endif
1723           if(me.eq.king.or..not.out1file)
1724      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1725
1726           if (isnan(etot) .or. etot.gt.1.0d4) then
1727             write (iout,*) "Energy too large",etot,
1728      &        " trying another random conformation"
1729             do itrial=1,100
1730               itmp=1
1731               call gen_rand_conf(itmp,*30)
1732               goto 40
1733    30         write (iout,*) 'Failed to generate random conformation',
1734      &          ', itrial=',itrial
1735               write (*,*) 'Processor:',me,
1736      &          ' Failed to generate random conformation',
1737      &          ' itrial=',itrial
1738               call intout
1739 #ifdef AIX
1740               call flush_(iout)
1741 #else
1742               call flush(iout)
1743 #endif
1744             enddo
1745             write (iout,'(a,i3,a)') 'Processor:',me,
1746      &        ' error in generating random conformation.'
1747             write (*,'(a,i3,a)') 'Processor:',me,
1748      &        ' error in generating random conformation.'
1749             call flush(iout)
1750 #ifdef MPI
1751             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1752 #else
1753             stop
1754 #endif
1755    40       continue
1756           else
1757             goto 44
1758           endif
1759           ENDDO
1760
1761           write (iout,'(a,i3,a)') 'Processor:',me,
1762      &        ' failed to generate a low-energy random conformation.'
1763             write (*,'(a,i3,a)') 'Processor:',me,
1764      &        ' failed to generate a low-energy random conformation.'
1765             call flush(iout)
1766 #ifdef MPI
1767             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1768 #else
1769             stop
1770 #endif
1771    44     continue
1772          endif
1773       endif
1774       call chainbuild_cart
1775       call kinetic(EK)
1776       if (tbf) then
1777         call verlet_bath
1778       endif
1779       kinetic_T=2.0d0/(dimen3*Rb)*EK
1780       if(me.eq.king.or..not.out1file)then
1781        call cartprint
1782        call intout
1783       endif
1784 #ifdef MPI
1785       tt0=MPI_Wtime()
1786 #else
1787       tt0=tcpu()
1788 #endif
1789       call zerograd
1790       call etotal(potEcomp)
1791       if (large) call enerprint(potEcomp)
1792 #ifdef TIMING_ENE
1793 #ifdef MPI
1794       t_etotal=t_etotal+MPI_Wtime()-tt0
1795 #else
1796       t_etotal=t_etotal+tcpu()-tt0
1797 #endif
1798 #endif
1799       potE=potEcomp(0)
1800       call cartgrad
1801       call lagrangian
1802       call max_accel
1803       if (amax*d_time .gt. dvmax) then
1804         d_time=d_time*dvmax/amax
1805         if(me.eq.king.or..not.out1file) write (iout,*)
1806      &   "Time step reduced to",d_time,
1807      &   " because of too large initial acceleration."
1808       endif
1809       if(me.eq.king.or..not.out1file)then
1810        write(iout,*) "Potential energy and its components"
1811        call enerprint(potEcomp)
1812 c       write(iout,*) (potEcomp(i),i=0,n_ene)
1813       endif
1814       potE=potEcomp(0)-potEcomp(20)
1815       totE=EK+potE
1816       itime=0
1817       if (ntwe.ne.0) call statout(itime)
1818       if(me.eq.king.or..not.out1file)
1819      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1820      &   " Kinetic energy",EK," Potential energy",potE,
1821      &   " Total energy",totE," Maximum acceleration ",
1822      &   amax
1823       if (large) then
1824         write (iout,*) "Initial coordinates"
1825         do i=1,nres
1826           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1827      &    (c(j,i+nres),j=1,3)
1828         enddo
1829         write (iout,*) "Initial dC"
1830         do i=0,nres
1831           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1832      &    (dc(j,i+nres),j=1,3)
1833         enddo
1834         write (iout,*) "Initial velocities"
1835         write (iout,"(13x,' backbone ',23x,' side chain')")
1836         do i=0,nres
1837           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1838      &    (d_t(j,i+nres),j=1,3)
1839         enddo
1840         write (iout,*) "Initial accelerations"
1841         do i=0,nres
1842 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1843           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1844      &    (d_a(j,i+nres),j=1,3)
1845         enddo
1846       endif
1847       do i=0,2*nres
1848         do j=1,3
1849           dc_old(j,i)=dc(j,i)
1850           d_t_old(j,i)=d_t(j,i)
1851           d_a_old(j,i)=d_a(j,i)
1852         enddo
1853 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1854       enddo
1855       if (RESPA) then
1856 #ifdef MPI
1857         tt0 =MPI_Wtime()
1858 #else
1859         tt0 = tcpu()
1860 #endif
1861         call zerograd
1862         call etotal_short(energia_short)
1863         if (large) call enerprint(potEcomp)
1864 #ifdef TIMING_ENE
1865 #ifdef MPI
1866         t_eshort=t_eshort+MPI_Wtime()-tt0
1867 #else
1868         t_eshort=t_eshort+tcpu()-tt0
1869 #endif
1870 #endif
1871         call cartgrad
1872         call lagrangian
1873         if(.not.out1file .and. large) then
1874           write (iout,*) "energia_long",energia_long(0),
1875      &     " energia_short",energia_short(0),
1876      &     " total",energia_long(0)+energia_short(0)
1877           write (iout,*) "Initial fast-force accelerations"
1878           do i=0,nres
1879             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1880      &      (d_a(j,i+nres),j=1,3)
1881           enddo
1882         endif
1883 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1884         do i=0,2*nres
1885           do j=1,3
1886             d_a_short(j,i)=d_a(j,i)
1887           enddo
1888         enddo
1889 #ifdef MPI
1890         tt0=MPI_Wtime()
1891 #else
1892         tt0=tcpu()
1893 #endif
1894         call zerograd
1895         call etotal_long(energia_long)
1896         if (large) call enerprint(potEcomp)
1897 #ifdef TIMING_ENE
1898 #ifdef MPI
1899         t_elong=t_elong+MPI_Wtime()-tt0
1900 #else
1901         t_elong=t_elong+tcpu()-tt0
1902 #endif
1903 #endif
1904         call cartgrad
1905         call lagrangian
1906         if(.not.out1file .and. large) then
1907           write (iout,*) "energia_long",energia_long(0)
1908           write (iout,*) "Initial slow-force accelerations"
1909           do i=0,nres
1910             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1911      &      (d_a(j,i+nres),j=1,3)
1912           enddo
1913         endif
1914 #ifdef MPI
1915         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1916 #else
1917         t_enegrad=t_enegrad+tcpu()-tt0
1918 #endif
1919       endif
1920       return
1921       end
1922 c-----------------------------------------------------------
1923       subroutine random_vel
1924       implicit real*8 (a-h,o-z)
1925       include 'DIMENSIONS'
1926       include 'COMMON.CONTROL'
1927       include 'COMMON.VAR'
1928       include 'COMMON.MD'
1929 #ifndef LANG0
1930       include 'COMMON.LANGEVIN'
1931 #else
1932       include 'COMMON.LANGEVIN.lang0'
1933 #endif
1934       include 'COMMON.CHAIN'
1935       include 'COMMON.DERIV'
1936       include 'COMMON.GEO'
1937       include 'COMMON.LOCAL'
1938       include 'COMMON.INTERACT'
1939       include 'COMMON.IOUNITS'
1940       include 'COMMON.NAMES'
1941       include 'COMMON.TIME1'
1942       double precision xv,sigv,lowb,highb,vec_afm(3)
1943 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1944 c First generate velocities in the eigenspace of the G matrix
1945 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1946 c      call flush(iout)
1947       xv=0.0d0
1948       ii=0
1949       do i=1,dimen
1950         do k=1,3
1951           ii=ii+1
1952           sigv=dsqrt((Rb*t_bath)/geigen(i))
1953           lowb=-5*sigv
1954           highb=5*sigv
1955           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1956
1957 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1958 c     &      " d_t_work_new",d_t_work_new(ii)
1959         enddo
1960       enddo
1961 C       if (SELFGUIDE.gt.0) then
1962 C       distance=0.0
1963 C       do j=1,3
1964 C       vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1965 C       distance=distance+vec_afm(j)**2
1966 C       enddo
1967 C       distance=dsqrt(distance)
1968 C       do j=1,3
1969 C         d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1970 C         d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1971 C         write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1972 C     &    d_t_work_new(j+(afmend-1)*3)
1973 C       enddo
1974
1975 C       endif
1976
1977 c diagnostics
1978 c      Ek1=0.0d0
1979 c      ii=0
1980 c      do i=1,dimen
1981 c        do k=1,3
1982 c          ii=ii+1
1983 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1984 c        enddo
1985 c      enddo
1986 c      write (iout,*) "Ek from eigenvectors",Ek1
1987 c end diagnostics
1988 c Transform velocities to UNRES coordinate space
1989       do k=0,2
1990         do i=1,dimen
1991           ind=(i-1)*3+k+1
1992           d_t_work(ind)=0.0d0
1993           do j=1,dimen
1994             d_t_work(ind)=d_t_work(ind)
1995      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1996           enddo
1997 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1998 c          call flush(iout)
1999         enddo
2000       enddo
2001 c Transfer to the d_t vector
2002       do j=1,3
2003         d_t(j,0)=d_t_work(j)
2004       enddo
2005       ind=3
2006       do i=nnt,nct-1
2007         do j=1,3
2008           ind=ind+1
2009           d_t(j,i)=d_t_work(ind)
2010         enddo
2011       enddo
2012       do i=nnt,nct
2013         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2014           do j=1,3
2015             ind=ind+1
2016             d_t(j,i+nres)=d_t_work(ind)
2017           enddo
2018         endif
2019       enddo
2020 c      call kinetic(EK)
2021 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2022 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2023 c      call flush(iout)
2024       return
2025       end
2026 #ifndef LANG0
2027 c-----------------------------------------------------------
2028       subroutine sd_verlet_p_setup
2029 c Sets up the parameters of stochastic Verlet algorithm
2030       implicit real*8 (a-h,o-z)
2031       include 'DIMENSIONS'
2032 #ifdef MPI
2033       include 'mpif.h'
2034 #endif
2035       include 'COMMON.CONTROL'
2036       include 'COMMON.VAR'
2037       include 'COMMON.MD'
2038 #ifndef LANG0
2039       include 'COMMON.LANGEVIN'
2040 #else
2041       include 'COMMON.LANGEVIN.lang0'
2042 #endif
2043       include 'COMMON.CHAIN'
2044       include 'COMMON.DERIV'
2045       include 'COMMON.GEO'
2046       include 'COMMON.LOCAL'
2047       include 'COMMON.INTERACT'
2048       include 'COMMON.IOUNITS'
2049       include 'COMMON.NAMES'
2050       include 'COMMON.TIME1'
2051       double precision emgdt(MAXRES6),
2052      & pterm,vterm,rho,rhoc,vsig,
2053      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2054      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2055      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2056       logical lprn /.false./
2057       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2058       double precision ktm
2059 #ifdef MPI
2060       tt0 = MPI_Wtime()
2061 #else
2062       tt0 = tcpu()
2063 #endif
2064 c
2065 c AL 8/17/04 Code adapted from tinker
2066 c
2067 c Get the frictional and random terms for stochastic dynamics in the
2068 c eigenspace of mass-scaled UNRES friction matrix
2069 c
2070       do i = 1, dimen
2071             gdt = fricgam(i) * d_time
2072 c
2073 c Stochastic dynamics reduces to simple MD for zero friction
2074 c
2075             if (gdt .le. zero) then
2076                pfric_vec(i) = 1.0d0
2077                vfric_vec(i) = d_time
2078                afric_vec(i) = 0.5d0 * d_time * d_time
2079                prand_vec(i) = 0.0d0
2080                vrand_vec1(i) = 0.0d0
2081                vrand_vec2(i) = 0.0d0
2082 c
2083 c Analytical expressions when friction coefficient is large
2084 c
2085             else
2086                if (gdt .ge. gdt_radius) then
2087                   egdt = dexp(-gdt)
2088                   pfric_vec(i) = egdt
2089                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2090                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2091                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2092                   vterm = 1.0d0 - egdt**2
2093                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2094 c
2095 c Use series expansions when friction coefficient is small
2096 c
2097                else
2098                   gdt2 = gdt * gdt
2099                   gdt3 = gdt * gdt2
2100                   gdt4 = gdt2 * gdt2
2101                   gdt5 = gdt2 * gdt3
2102                   gdt6 = gdt3 * gdt3
2103                   gdt7 = gdt3 * gdt4
2104                   gdt8 = gdt4 * gdt4
2105                   gdt9 = gdt4 * gdt5
2106                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2107      &                          - gdt5/120.0d0 + gdt6/720.0d0
2108      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2109      &                          - gdt9/362880.0d0) / fricgam(i)**2
2110                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2111                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2112                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2113      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2114      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2115      &                       + 127.0d0*gdt9/90720.0d0
2116                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2117      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2118      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2119      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2120                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2121      &                       - 17.0d0*gdt2/1280.0d0
2122      &                       + 17.0d0*gdt3/6144.0d0
2123      &                       + 40967.0d0*gdt4/34406400.0d0
2124      &                       - 57203.0d0*gdt5/275251200.0d0
2125      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2126                end if
2127 c
2128 c Compute the scaling factors of random terms for the nonzero friction case
2129 c
2130                ktm = 0.5d0*d_time/fricgam(i)
2131                psig = dsqrt(ktm*pterm) / fricgam(i)
2132                vsig = dsqrt(ktm*vterm)
2133                rhoc = dsqrt(1.0d0 - rho*rho)
2134                prand_vec(i) = psig
2135                vrand_vec1(i) = vsig * rho
2136                vrand_vec2(i) = vsig * rhoc
2137             end if
2138       end do
2139       if (lprn) then
2140       write (iout,*)
2141      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2142      &  " vrand_vec2"
2143       do i=1,dimen
2144         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2145      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2146       enddo
2147       endif
2148 c
2149 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2150 c
2151 #ifndef   LANG0
2152       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2153       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2154       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2155       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2156       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2157       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2158 #endif
2159 #ifdef MPI
2160       t_sdsetup=t_sdsetup+MPI_Wtime()
2161 #else
2162       t_sdsetup=t_sdsetup+tcpu()-tt0
2163 #endif
2164       return
2165       end
2166 c-------------------------------------------------------------
2167       subroutine eigtransf1(n,ndim,ab,d,c)
2168       implicit none
2169       integer n,ndim
2170       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2171       integer i,j,k
2172       do i=1,n
2173         do j=1,n
2174           c(i,j)=0.0d0
2175           do k=1,n
2176             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2177           enddo
2178         enddo
2179       enddo
2180       return
2181       end
2182 c-------------------------------------------------------------
2183       subroutine eigtransf(n,ndim,a,b,d,c)
2184       implicit none
2185       integer n,ndim
2186       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2187       integer i,j,k
2188       do i=1,n
2189         do j=1,n
2190           c(i,j)=0.0d0
2191           do k=1,n
2192             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2193           enddo
2194         enddo
2195       enddo
2196       return
2197       end
2198 c-------------------------------------------------------------
2199       subroutine sd_verlet1
2200 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2201       implicit real*8 (a-h,o-z)
2202       include 'DIMENSIONS'
2203       include 'COMMON.CONTROL'
2204       include 'COMMON.VAR'
2205       include 'COMMON.MD'
2206 #ifndef LANG0
2207       include 'COMMON.LANGEVIN'
2208 #else
2209       include 'COMMON.LANGEVIN.lang0'
2210 #endif
2211       include 'COMMON.CHAIN'
2212       include 'COMMON.DERIV'
2213       include 'COMMON.GEO'
2214       include 'COMMON.LOCAL'
2215       include 'COMMON.INTERACT'
2216       include 'COMMON.IOUNITS'
2217       include 'COMMON.NAMES'
2218       double precision stochforcvec(MAXRES6)
2219       common /stochcalc/ stochforcvec
2220       logical lprn /.false./
2221
2222 c      write (iout,*) "dc_old"
2223 c      do i=0,nres
2224 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2225 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2226 c      enddo
2227       do j=1,3
2228         dc_work(j)=dc_old(j,0)
2229         d_t_work(j)=d_t_old(j,0)
2230         d_a_work(j)=d_a_old(j,0)
2231       enddo
2232       ind=3
2233       do i=nnt,nct-1
2234         do j=1,3
2235           dc_work(ind+j)=dc_old(j,i)
2236           d_t_work(ind+j)=d_t_old(j,i)
2237           d_a_work(ind+j)=d_a_old(j,i)
2238         enddo
2239         ind=ind+3
2240       enddo
2241       do i=nnt,nct
2242         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2243           do j=1,3
2244             dc_work(ind+j)=dc_old(j,i+nres)
2245             d_t_work(ind+j)=d_t_old(j,i+nres)
2246             d_a_work(ind+j)=d_a_old(j,i+nres)
2247           enddo
2248           ind=ind+3
2249         endif
2250       enddo
2251 #ifndef LANG0
2252       if (lprn) then
2253       write (iout,*)
2254      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2255      &  " vrand_mat2"
2256       do i=1,dimen
2257         do j=1,dimen
2258           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2259      &      vfric_mat(i,j),afric_mat(i,j),
2260      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2261         enddo
2262       enddo
2263       endif
2264       do i=1,dimen
2265         ddt1=0.0d0
2266         ddt2=0.0d0
2267         do j=1,dimen
2268           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2269      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2270           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2271           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2272         enddo
2273         d_t_work_new(i)=ddt1+0.5d0*ddt2
2274         d_t_work(i)=ddt1+ddt2
2275       enddo
2276 #endif
2277       do j=1,3
2278         dc(j,0)=dc_work(j)
2279         d_t(j,0)=d_t_work(j)
2280       enddo
2281       ind=3
2282       do i=nnt,nct-1
2283         do j=1,3
2284           dc(j,i)=dc_work(ind+j)
2285           d_t(j,i)=d_t_work(ind+j)
2286         enddo
2287         ind=ind+3
2288       enddo
2289       do i=nnt,nct
2290         if (itype(i).ne.10) then
2291           inres=i+nres
2292           do j=1,3
2293             dc(j,inres)=dc_work(ind+j)
2294             d_t(j,inres)=d_t_work(ind+j)
2295           enddo
2296           ind=ind+3
2297         endif
2298       enddo
2299       return
2300       end
2301 c--------------------------------------------------------------------------
2302       subroutine sd_verlet2
2303 c  Calculating the adjusted velocities for accelerations
2304       implicit real*8 (a-h,o-z)
2305       include 'DIMENSIONS'
2306       include 'COMMON.CONTROL'
2307       include 'COMMON.VAR'
2308       include 'COMMON.MD'
2309 #ifndef LANG0
2310       include 'COMMON.LANGEVIN'
2311 #else
2312       include 'COMMON.LANGEVIN.lang0'
2313 #endif
2314       include 'COMMON.CHAIN'
2315       include 'COMMON.DERIV'
2316       include 'COMMON.GEO'
2317       include 'COMMON.LOCAL'
2318       include 'COMMON.INTERACT'
2319       include 'COMMON.IOUNITS'
2320       include 'COMMON.NAMES'
2321       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2322       common /stochcalc/ stochforcvec
2323 c
2324 c Compute the stochastic forces which contribute to velocity change
2325 c
2326       call stochastic_force(stochforcvecV)
2327
2328 #ifndef LANG0
2329       do i=1,dimen
2330         ddt1=0.0d0
2331         ddt2=0.0d0
2332         do j=1,dimen
2333           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2334           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2335      &     vrand_mat2(i,j)*stochforcvecV(j)
2336         enddo
2337         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2338       enddo
2339 #endif
2340       do j=1,3
2341         d_t(j,0)=d_t_work(j)
2342       enddo
2343       ind=3
2344       do i=nnt,nct-1
2345         do j=1,3
2346           d_t(j,i)=d_t_work(ind+j)
2347         enddo
2348         ind=ind+3
2349       enddo
2350       do i=nnt,nct
2351         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2352           inres=i+nres
2353           do j=1,3
2354             d_t(j,inres)=d_t_work(ind+j)
2355           enddo
2356           ind=ind+3
2357         endif
2358       enddo
2359       return
2360       end
2361 c-----------------------------------------------------------
2362       subroutine sd_verlet_ciccotti_setup
2363 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2364 c version
2365       implicit real*8 (a-h,o-z)
2366       include 'DIMENSIONS'
2367 #ifdef MPI
2368       include 'mpif.h'
2369 #endif
2370       include 'COMMON.CONTROL'
2371       include 'COMMON.VAR'
2372       include 'COMMON.MD'
2373 #ifndef LANG0
2374       include 'COMMON.LANGEVIN'
2375 #else
2376       include 'COMMON.LANGEVIN.lang0'
2377 #endif
2378       include 'COMMON.CHAIN'
2379       include 'COMMON.DERIV'
2380       include 'COMMON.GEO'
2381       include 'COMMON.LOCAL'
2382       include 'COMMON.INTERACT'
2383       include 'COMMON.IOUNITS'
2384       include 'COMMON.NAMES'
2385       include 'COMMON.TIME1'
2386       double precision emgdt(MAXRES6),
2387      & pterm,vterm,rho,rhoc,vsig,
2388      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2389      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2390      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2391       logical lprn /.false./
2392       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2393       double precision ktm
2394 #ifdef MPI
2395       tt0 = MPI_Wtime()
2396 #else
2397       tt0 = tcpu()
2398 #endif
2399 c
2400 c AL 8/17/04 Code adapted from tinker
2401 c
2402 c Get the frictional and random terms for stochastic dynamics in the
2403 c eigenspace of mass-scaled UNRES friction matrix
2404 c
2405       do i = 1, dimen
2406             write (iout,*) "i",i," fricgam",fricgam(i)
2407             gdt = fricgam(i) * d_time
2408 c
2409 c Stochastic dynamics reduces to simple MD for zero friction
2410 c
2411             if (gdt .le. zero) then
2412                pfric_vec(i) = 1.0d0
2413                vfric_vec(i) = d_time
2414                afric_vec(i) = 0.5d0*d_time*d_time
2415                prand_vec(i) = afric_vec(i)
2416                vrand_vec2(i) = vfric_vec(i)
2417 c
2418 c Analytical expressions when friction coefficient is large
2419 c
2420             else
2421                egdt = dexp(-gdt)
2422                pfric_vec(i) = egdt
2423                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2424                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2425                prand_vec(i) = afric_vec(i)
2426                vrand_vec2(i) = vfric_vec(i)
2427 c
2428 c Compute the scaling factors of random terms for the nonzero friction case
2429 c
2430 c               ktm = 0.5d0*d_time/fricgam(i)
2431 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2432 c               vsig = dsqrt(ktm*vterm)
2433 c               prand_vec(i) = psig*afric_vec(i)
2434 c               vrand_vec2(i) = vsig*vfric_vec(i)
2435             end if
2436       end do
2437       if (lprn) then
2438       write (iout,*)
2439      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2440      &  " vrand_vec2"
2441       do i=1,dimen
2442         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2443      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2444       enddo
2445       endif
2446 c
2447 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2448 c
2449       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2450       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2451       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2452       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2453       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2454 #ifdef MPI
2455       t_sdsetup=t_sdsetup+MPI_Wtime()
2456 #else
2457       t_sdsetup=t_sdsetup+tcpu()-tt0
2458 #endif
2459       return
2460       end
2461 c-------------------------------------------------------------
2462       subroutine sd_verlet1_ciccotti
2463 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2464       implicit real*8 (a-h,o-z)
2465       include 'DIMENSIONS'
2466 #ifdef MPI
2467       include 'mpif.h'
2468 #endif
2469       include 'COMMON.CONTROL'
2470       include 'COMMON.VAR'
2471       include 'COMMON.MD'
2472 #ifndef LANG0
2473       include 'COMMON.LANGEVIN'
2474 #else
2475       include 'COMMON.LANGEVIN.lang0'
2476 #endif
2477       include 'COMMON.CHAIN'
2478       include 'COMMON.DERIV'
2479       include 'COMMON.GEO'
2480       include 'COMMON.LOCAL'
2481       include 'COMMON.INTERACT'
2482       include 'COMMON.IOUNITS'
2483       include 'COMMON.NAMES'
2484       double precision stochforcvec(MAXRES6)
2485       common /stochcalc/ stochforcvec
2486       logical lprn /.false./
2487
2488 c      write (iout,*) "dc_old"
2489 c      do i=0,nres
2490 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2491 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2492 c      enddo
2493       do j=1,3
2494         dc_work(j)=dc_old(j,0)
2495         d_t_work(j)=d_t_old(j,0)
2496         d_a_work(j)=d_a_old(j,0)
2497       enddo
2498       ind=3
2499       do i=nnt,nct-1
2500         do j=1,3
2501           dc_work(ind+j)=dc_old(j,i)
2502           d_t_work(ind+j)=d_t_old(j,i)
2503           d_a_work(ind+j)=d_a_old(j,i)
2504         enddo
2505         ind=ind+3
2506       enddo
2507       do i=nnt,nct
2508         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2509           do j=1,3
2510             dc_work(ind+j)=dc_old(j,i+nres)
2511             d_t_work(ind+j)=d_t_old(j,i+nres)
2512             d_a_work(ind+j)=d_a_old(j,i+nres)
2513           enddo
2514           ind=ind+3
2515         endif
2516       enddo
2517
2518 #ifndef LANG0
2519       if (lprn) then
2520       write (iout,*)
2521      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2522      &  " vrand_mat2"
2523       do i=1,dimen
2524         do j=1,dimen
2525           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2526      &      vfric_mat(i,j),afric_mat(i,j),
2527      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2528         enddo
2529       enddo
2530       endif
2531       do i=1,dimen
2532         ddt1=0.0d0
2533         ddt2=0.0d0
2534         do j=1,dimen
2535           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2536      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2537           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2538           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2539         enddo
2540         d_t_work_new(i)=ddt1+0.5d0*ddt2
2541         d_t_work(i)=ddt1+ddt2
2542       enddo
2543 #endif
2544       do j=1,3
2545         dc(j,0)=dc_work(j)
2546         d_t(j,0)=d_t_work(j)
2547       enddo
2548       ind=3
2549       do i=nnt,nct-1
2550         do j=1,3
2551           dc(j,i)=dc_work(ind+j)
2552           d_t(j,i)=d_t_work(ind+j)
2553         enddo
2554         ind=ind+3
2555       enddo
2556       do i=nnt,nct
2557         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2558           inres=i+nres
2559           do j=1,3
2560             dc(j,inres)=dc_work(ind+j)
2561             d_t(j,inres)=d_t_work(ind+j)
2562           enddo
2563           ind=ind+3
2564         endif
2565       enddo
2566       return
2567       end
2568 c--------------------------------------------------------------------------
2569       subroutine sd_verlet2_ciccotti
2570 c  Calculating the adjusted velocities for accelerations
2571       implicit real*8 (a-h,o-z)
2572       include 'DIMENSIONS'
2573       include 'COMMON.CONTROL'
2574       include 'COMMON.VAR'
2575       include 'COMMON.MD'
2576 #ifndef LANG0
2577       include 'COMMON.LANGEVIN'
2578 #else
2579       include 'COMMON.LANGEVIN.lang0'
2580 #endif
2581       include 'COMMON.CHAIN'
2582       include 'COMMON.DERIV'
2583       include 'COMMON.GEO'
2584       include 'COMMON.LOCAL'
2585       include 'COMMON.INTERACT'
2586       include 'COMMON.IOUNITS'
2587       include 'COMMON.NAMES'
2588       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2589       common /stochcalc/ stochforcvec
2590 c
2591 c Compute the stochastic forces which contribute to velocity change
2592 c
2593       call stochastic_force(stochforcvecV)
2594 #ifndef LANG0
2595       do i=1,dimen
2596         ddt1=0.0d0
2597         ddt2=0.0d0
2598         do j=1,dimen
2599
2600           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2601 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2602           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2603         enddo
2604         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2605       enddo
2606 #endif
2607       do j=1,3
2608         d_t(j,0)=d_t_work(j)
2609       enddo
2610       ind=3
2611       do i=nnt,nct-1
2612         do j=1,3
2613           d_t(j,i)=d_t_work(ind+j)
2614         enddo
2615         ind=ind+3
2616       enddo
2617       do i=nnt,nct
2618         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2619           inres=i+nres
2620           do j=1,3
2621             d_t(j,inres)=d_t_work(ind+j)
2622           enddo
2623           ind=ind+3
2624         endif
2625       enddo
2626       return
2627       end
2628 #endif