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           if (itype(i).ne.ntyp1) 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 #if defined(AIX) || defined(PGI)
1503         open(ipdb,
1504      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1505      &  //".pdb", position='append')
1506 #else
1507         open(ipdb,
1508      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1509      &  //".pdb", access='append')
1510 #endif
1511       else
1512 #ifdef NOXDR
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))//".x")
1516         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1517      &  //".x"
1518 #else
1519         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1520      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1521      &      liczba(:ilen(liczba))//".cx")
1522         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1523      &  //".cx"
1524 #endif
1525       endif
1526 #else
1527       if(mdpdb) then
1528          if (ilen(tmpdir).gt.0)
1529      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1530          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1531       else
1532          if (ilen(tmpdir).gt.0)
1533      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1534          cartname=prefix(:ilen(prefix))//"_MD.cx"
1535       endif
1536 #endif
1537       if (usampl) then
1538         write (qstr,'(256(1h ))')
1539         ipos=1
1540         do i=1,nfrag
1541           iq = qinfrag(i,iset)*10
1542           iw = wfrag(i,iset)/100
1543           if (iw.gt.0) then
1544             if(me.eq.king.or..not.out1file)
1545      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1546             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1547             ipos=ipos+7
1548           endif
1549         enddo
1550         do i=1,npair
1551           iq = qinpair(i,iset)*10
1552           iw = wpair(i,iset)/100
1553           if (iw.gt.0) then
1554             if(me.eq.king.or..not.out1file)
1555      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1556             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1557             ipos=ipos+7
1558           endif
1559         enddo
1560 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1561 #ifdef NOXDR
1562 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1563 #else
1564 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1565 #endif
1566 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1567       endif
1568       icg=1
1569       if (rest) then
1570        if (restart1file) then
1571          if (me.eq.king)
1572      &     inquire(file=mremd_rst_name,exist=file_exist)
1573 #ifdef MPI
1574            write (*,*) me," Before broadcast: file_exist",file_exist
1575          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1576      &          IERR)
1577          write (*,*) me," After broadcast: file_exist",file_exist
1578 c        inquire(file=mremd_rst_name,exist=file_exist)
1579 #endif
1580         if(me.eq.king.or..not.out1file)
1581      &   write(iout,*) "Initial state read by master and distributed"
1582        else
1583          if (ilen(tmpdir).gt.0)
1584      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1585      &      //liczba(:ilen(liczba))//'.rst')
1586         inquire(file=rest2name,exist=file_exist)
1587        endif
1588        if(file_exist) then
1589          if(.not.restart1file) then
1590            if(me.eq.king.or..not.out1file)
1591      &      write(iout,*) "Initial state will be read from file ",
1592      &      rest2name(:ilen(rest2name))
1593            call readrst
1594          endif
1595          call rescale_weights(t_bath)
1596        else
1597         if(me.eq.king.or..not.out1file)then
1598          if (restart1file) then
1599           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1600      &       " does not exist"
1601          else
1602           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1603      &       " does not exist"
1604          endif
1605          write(iout,*) "Initial velocities randomly generated"
1606         endif
1607         call random_vel
1608         totT=0.0d0
1609         totTafm=totT
1610        endif
1611       else
1612 c Generate initial velocities
1613         if(me.eq.king.or..not.out1file)
1614      &   write(iout,*) "Initial velocities randomly generated"
1615         call random_vel
1616         totT=0.0d0
1617 CtotTafm is the variable for AFM time which eclipsed during
1618         totTafm=totT
1619       endif
1620 c      rest2name = prefix(:ilen(prefix))//'.rst'
1621       if(me.eq.king.or..not.out1file)then
1622        write (iout,*) "Initial velocities"
1623        do i=0,nres
1624          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1625      &   (d_t(j,i+nres),j=1,3)
1626        enddo
1627 c  Zeroing the total angular momentum of the system
1628        write(iout,*) "Calling the zero-angular
1629      & momentum subroutine"
1630       endif
1631       call inertia_tensor
1632 c  Getting the potential energy and forces and velocities and accelerations
1633       call vcm_vel(vcm)
1634 c      write (iout,*) "velocity of the center of the mass:"
1635 c      write (iout,*) (vcm(j),j=1,3)
1636       do j=1,3
1637         d_t(j,0)=d_t(j,0)-vcm(j)
1638       enddo
1639 c Removing the velocity of the center of mass
1640       call vcm_vel(vcm)
1641       if(me.eq.king.or..not.out1file)then
1642        write (iout,*) "vcm right after adjustment:"
1643        write (iout,*) (vcm(j),j=1,3)
1644       endif
1645       if (.not.rest) then
1646          if (start_from_model) then
1647           i_model=iran_num(1,constr_homology)
1648           write (iout,*) 'starting from model ',i_model
1649           do i=1,2*nres
1650            do j=1,3
1651             c(j,i)=chomo(j,i,i_model)
1652            enddo
1653           enddo
1654
1655           call int_from_cart(.true.,.false.)
1656           call sc_loc_geom(.false.)
1657           do i=1,nres-1
1658             do j=1,3
1659              dc(j,i)=c(j,i+1)-c(j,i)
1660              dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1661             enddo
1662           enddo
1663           do i=2,nres-1
1664             do j=1,3
1665              dc(j,i+nres)=c(j,i+nres)-c(j,i)
1666              dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1667             enddo
1668           enddo
1669           call etotal(energia(0))
1670          endif
1671 c         call chainbuild_cart
1672          write (iout,*) "Initial energies"
1673          call enerprint(energia(0))
1674          write (iout,*) "PREMINIM ",preminim
1675          if(iranconf.ne.0 .or. preminim) then
1676           if (overlapsc) then
1677            write (iout,*) 'Calling OVERLAP_SC'
1678            call overlap_sc(fail)
1679           endif
1680
1681           if (searchsc) then
1682            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1683            print *,'SC_move',nft_sc,etot
1684            if(me.eq.king.or..not.out1file)
1685      &      write(iout,*) 'SC_move',nft_sc,etot
1686           endif
1687
1688           if(dccart)then
1689            print *, 'Calling MINIM_DC'
1690            call minim_dc(etot,iretcode,nfun)
1691           else
1692            call geom_to_var(nvar,varia)
1693            print *,'Calling MINIMIZE.'
1694            call minimize(etot,varia,iretcode,nfun)
1695            call var_to_geom(nvar,varia)
1696           endif
1697           if(me.eq.king.or..not.out1file) then
1698              write(iout,*) "Minimized energy is",etot
1699              write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1700              call etotal(potEcomp)
1701              call enerprint(potEcomp)
1702           endif
1703          endif
1704          if(iranconf.ne.0) then
1705 c 8/22/17 AL Loop to produce a low-energy random conformation
1706           DO iranmin=1,10
1707           if (overlapsc) then
1708            if(me.eq.king.or..not.out1file)
1709      &        write (iout,*) 'Calling OVERLAP_SC'
1710            call overlap_sc(fail)
1711           endif
1712
1713           if (searchsc) then
1714            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1715            print *,'SC_move',nft_sc,etot
1716            if(me.eq.king.or..not.out1file)
1717      &      write(iout,*) 'SC_move',nft_sc,etot
1718           endif
1719
1720           if(dccart)then
1721            print *, 'Calling MINIM_DC'
1722            call minim_dc(etot,iretcode,nfun)
1723           else
1724            call geom_to_var(nvar,varia)
1725            print *,'Calling MINIMIZE.'
1726            call minimize(etot,varia,iretcode,nfun)
1727            call var_to_geom(nvar,varia)
1728           endif
1729           if(me.eq.king.or..not.out1file)
1730      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1731
1732           if (isnan(etot) .or. etot.gt.1.0d4) then
1733             write (iout,*) "Energy too large",etot,
1734      &        " trying another random conformation"
1735             do itrial=1,100
1736               itmp=1
1737               call gen_rand_conf(itmp,*30)
1738               goto 40
1739    30         write (iout,*) 'Failed to generate random conformation',
1740      &          ', itrial=',itrial
1741               write (*,*) 'Processor:',me,
1742      &          ' Failed to generate random conformation',
1743      &          ' itrial=',itrial
1744               call intout
1745 #ifdef AIX
1746               call flush_(iout)
1747 #else
1748               call flush(iout)
1749 #endif
1750             enddo
1751             write (iout,'(a,i3,a)') 'Processor:',me,
1752      &        ' error in generating random conformation.'
1753             write (*,'(a,i3,a)') 'Processor:',me,
1754      &        ' error in generating random conformation.'
1755             call flush(iout)
1756 #ifdef MPI
1757             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1758 #else
1759             stop
1760 #endif
1761    40       continue
1762           else
1763             goto 44
1764           endif
1765           ENDDO
1766
1767           write (iout,'(a,i3,a)') 'Processor:',me,
1768      &        ' failed to generate a low-energy random conformation.'
1769             write (*,'(a,i3,a)') 'Processor:',me,
1770      &        ' failed to generate a low-energy random conformation.'
1771             call flush(iout)
1772 #ifdef MPI
1773             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1774 #else
1775             stop
1776 #endif
1777    44     continue
1778          endif
1779       endif
1780       call chainbuild_cart
1781       call kinetic(EK)
1782       if (tbf) then
1783         call verlet_bath
1784       endif
1785       kinetic_T=2.0d0/(dimen3*Rb)*EK
1786       if(me.eq.king.or..not.out1file)then
1787        call cartprint
1788        call intout
1789       endif
1790 #ifdef MPI
1791       tt0=MPI_Wtime()
1792 #else
1793       tt0=tcpu()
1794 #endif
1795       call zerograd
1796       call etotal(potEcomp)
1797       if (large) call enerprint(potEcomp)
1798 #ifdef TIMING_ENE
1799 #ifdef MPI
1800       t_etotal=t_etotal+MPI_Wtime()-tt0
1801 #else
1802       t_etotal=t_etotal+tcpu()-tt0
1803 #endif
1804 #endif
1805       potE=potEcomp(0)
1806       call cartgrad
1807       call lagrangian
1808       call max_accel
1809       if (amax*d_time .gt. dvmax) then
1810         d_time=d_time*dvmax/amax
1811         if(me.eq.king.or..not.out1file) write (iout,*)
1812      &   "Time step reduced to",d_time,
1813      &   " because of too large initial acceleration."
1814       endif
1815       if(me.eq.king.or..not.out1file)then
1816        write(iout,*) "Potential energy and its components"
1817        call enerprint(potEcomp)
1818 c       write(iout,*) (potEcomp(i),i=0,n_ene)
1819       endif
1820       potE=potEcomp(0)-potEcomp(20)
1821       totE=EK+potE
1822       itime=0
1823       if (ntwe.ne.0) call statout(itime)
1824       if(me.eq.king.or..not.out1file)
1825      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1826      &   " Kinetic energy",EK," Potential energy",potE,
1827      &   " Total energy",totE," Maximum acceleration ",
1828      &   amax
1829       if (large) then
1830         write (iout,*) "Initial coordinates"
1831         do i=1,nres
1832           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1833      &    (c(j,i+nres),j=1,3)
1834         enddo
1835         write (iout,*) "Initial dC"
1836         do i=0,nres
1837           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1838      &    (dc(j,i+nres),j=1,3)
1839         enddo
1840         write (iout,*) "Initial velocities"
1841         write (iout,"(13x,' backbone ',23x,' side chain')")
1842         do i=0,nres
1843           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1844      &    (d_t(j,i+nres),j=1,3)
1845         enddo
1846         write (iout,*) "Initial accelerations"
1847         do i=0,nres
1848 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1849           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1850      &    (d_a(j,i+nres),j=1,3)
1851         enddo
1852       endif
1853       do i=0,2*nres
1854         do j=1,3
1855           dc_old(j,i)=dc(j,i)
1856           d_t_old(j,i)=d_t(j,i)
1857           d_a_old(j,i)=d_a(j,i)
1858         enddo
1859 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1860       enddo
1861       if (RESPA) then
1862 #ifdef MPI
1863         tt0 =MPI_Wtime()
1864 #else
1865         tt0 = tcpu()
1866 #endif
1867         call zerograd
1868         call etotal_short(energia_short)
1869         if (large) call enerprint(potEcomp)
1870 #ifdef TIMING_ENE
1871 #ifdef MPI
1872         t_eshort=t_eshort+MPI_Wtime()-tt0
1873 #else
1874         t_eshort=t_eshort+tcpu()-tt0
1875 #endif
1876 #endif
1877         call cartgrad
1878         call lagrangian
1879         if(.not.out1file .and. large) then
1880           write (iout,*) "energia_long",energia_long(0),
1881      &     " energia_short",energia_short(0),
1882      &     " total",energia_long(0)+energia_short(0)
1883           write (iout,*) "Initial fast-force accelerations"
1884           do i=0,nres
1885             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1886      &      (d_a(j,i+nres),j=1,3)
1887           enddo
1888         endif
1889 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1890         do i=0,2*nres
1891           do j=1,3
1892             d_a_short(j,i)=d_a(j,i)
1893           enddo
1894         enddo
1895 #ifdef MPI
1896         tt0=MPI_Wtime()
1897 #else
1898         tt0=tcpu()
1899 #endif
1900         call zerograd
1901         call etotal_long(energia_long)
1902         if (large) call enerprint(potEcomp)
1903 #ifdef TIMING_ENE
1904 #ifdef MPI
1905         t_elong=t_elong+MPI_Wtime()-tt0
1906 #else
1907         t_elong=t_elong+tcpu()-tt0
1908 #endif
1909 #endif
1910         call cartgrad
1911         call lagrangian
1912         if(.not.out1file .and. large) then
1913           write (iout,*) "energia_long",energia_long(0)
1914           write (iout,*) "Initial slow-force accelerations"
1915           do i=0,nres
1916             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1917      &      (d_a(j,i+nres),j=1,3)
1918           enddo
1919         endif
1920 #ifdef MPI
1921         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1922 #else
1923         t_enegrad=t_enegrad+tcpu()-tt0
1924 #endif
1925       endif
1926       return
1927       end
1928 c-----------------------------------------------------------
1929       subroutine random_vel
1930       implicit real*8 (a-h,o-z)
1931       include 'DIMENSIONS'
1932       include 'COMMON.CONTROL'
1933       include 'COMMON.VAR'
1934       include 'COMMON.MD'
1935 #ifndef LANG0
1936       include 'COMMON.LANGEVIN'
1937 #else
1938       include 'COMMON.LANGEVIN.lang0'
1939 #endif
1940       include 'COMMON.CHAIN'
1941       include 'COMMON.DERIV'
1942       include 'COMMON.GEO'
1943       include 'COMMON.LOCAL'
1944       include 'COMMON.INTERACT'
1945       include 'COMMON.IOUNITS'
1946       include 'COMMON.NAMES'
1947       include 'COMMON.TIME1'
1948       double precision xv,sigv,lowb,highb,vec_afm(3)
1949 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1950 c First generate velocities in the eigenspace of the G matrix
1951 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1952 c      call flush(iout)
1953       xv=0.0d0
1954       ii=0
1955       do i=1,dimen
1956         do k=1,3
1957           ii=ii+1
1958           sigv=dsqrt((Rb*t_bath)/geigen(i))
1959           lowb=-5*sigv
1960           highb=5*sigv
1961           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1962
1963 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1964 c     &      " d_t_work_new",d_t_work_new(ii)
1965         enddo
1966       enddo
1967 C       if (SELFGUIDE.gt.0) then
1968 C       distance=0.0
1969 C       do j=1,3
1970 C       vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1971 C       distance=distance+vec_afm(j)**2
1972 C       enddo
1973 C       distance=dsqrt(distance)
1974 C       do j=1,3
1975 C         d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1976 C         d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1977 C         write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1978 C     &    d_t_work_new(j+(afmend-1)*3)
1979 C       enddo
1980
1981 C       endif
1982
1983 c diagnostics
1984 c      Ek1=0.0d0
1985 c      ii=0
1986 c      do i=1,dimen
1987 c        do k=1,3
1988 c          ii=ii+1
1989 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1990 c        enddo
1991 c      enddo
1992 c      write (iout,*) "Ek from eigenvectors",Ek1
1993 c end diagnostics
1994 c Transform velocities to UNRES coordinate space
1995       do k=0,2
1996         do i=1,dimen
1997           ind=(i-1)*3+k+1
1998           d_t_work(ind)=0.0d0
1999           do j=1,dimen
2000             d_t_work(ind)=d_t_work(ind)
2001      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2002           enddo
2003 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2004 c          call flush(iout)
2005         enddo
2006       enddo
2007 c Transfer to the d_t vector
2008       do j=1,3
2009         d_t(j,0)=d_t_work(j)
2010       enddo
2011       ind=3
2012       do i=nnt,nct-1
2013         do j=1,3
2014           ind=ind+1
2015           d_t(j,i)=d_t_work(ind)
2016         enddo
2017       enddo
2018       do i=nnt,nct
2019         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2020           do j=1,3
2021             ind=ind+1
2022             d_t(j,i+nres)=d_t_work(ind)
2023           enddo
2024         endif
2025       enddo
2026 c      call kinetic(EK)
2027 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2028 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2029 c      call flush(iout)
2030       return
2031       end
2032 #ifndef LANG0
2033 c-----------------------------------------------------------
2034       subroutine sd_verlet_p_setup
2035 c Sets up the parameters of stochastic Verlet algorithm
2036       implicit real*8 (a-h,o-z)
2037       include 'DIMENSIONS'
2038 #ifdef MPI
2039       include 'mpif.h'
2040 #endif
2041       include 'COMMON.CONTROL'
2042       include 'COMMON.VAR'
2043       include 'COMMON.MD'
2044 #ifndef LANG0
2045       include 'COMMON.LANGEVIN'
2046 #else
2047       include 'COMMON.LANGEVIN.lang0'
2048 #endif
2049       include 'COMMON.CHAIN'
2050       include 'COMMON.DERIV'
2051       include 'COMMON.GEO'
2052       include 'COMMON.LOCAL'
2053       include 'COMMON.INTERACT'
2054       include 'COMMON.IOUNITS'
2055       include 'COMMON.NAMES'
2056       include 'COMMON.TIME1'
2057       double precision emgdt(MAXRES6),
2058      & pterm,vterm,rho,rhoc,vsig,
2059      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2060      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2061      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2062       logical lprn /.false./
2063       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2064       double precision ktm
2065 #ifdef MPI
2066       tt0 = MPI_Wtime()
2067 #else
2068       tt0 = tcpu()
2069 #endif
2070 c
2071 c AL 8/17/04 Code adapted from tinker
2072 c
2073 c Get the frictional and random terms for stochastic dynamics in the
2074 c eigenspace of mass-scaled UNRES friction matrix
2075 c
2076       do i = 1, dimen
2077             gdt = fricgam(i) * d_time
2078 c
2079 c Stochastic dynamics reduces to simple MD for zero friction
2080 c
2081             if (gdt .le. zero) then
2082                pfric_vec(i) = 1.0d0
2083                vfric_vec(i) = d_time
2084                afric_vec(i) = 0.5d0 * d_time * d_time
2085                prand_vec(i) = 0.0d0
2086                vrand_vec1(i) = 0.0d0
2087                vrand_vec2(i) = 0.0d0
2088 c
2089 c Analytical expressions when friction coefficient is large
2090 c
2091             else
2092                if (gdt .ge. gdt_radius) then
2093                   egdt = dexp(-gdt)
2094                   pfric_vec(i) = egdt
2095                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2096                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2097                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2098                   vterm = 1.0d0 - egdt**2
2099                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2100 c
2101 c Use series expansions when friction coefficient is small
2102 c
2103                else
2104                   gdt2 = gdt * gdt
2105                   gdt3 = gdt * gdt2
2106                   gdt4 = gdt2 * gdt2
2107                   gdt5 = gdt2 * gdt3
2108                   gdt6 = gdt3 * gdt3
2109                   gdt7 = gdt3 * gdt4
2110                   gdt8 = gdt4 * gdt4
2111                   gdt9 = gdt4 * gdt5
2112                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2113      &                          - gdt5/120.0d0 + gdt6/720.0d0
2114      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2115      &                          - gdt9/362880.0d0) / fricgam(i)**2
2116                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2117                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2118                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2119      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2120      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2121      &                       + 127.0d0*gdt9/90720.0d0
2122                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2123      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2124      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2125      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2126                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2127      &                       - 17.0d0*gdt2/1280.0d0
2128      &                       + 17.0d0*gdt3/6144.0d0
2129      &                       + 40967.0d0*gdt4/34406400.0d0
2130      &                       - 57203.0d0*gdt5/275251200.0d0
2131      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2132                end if
2133 c
2134 c Compute the scaling factors of random terms for the nonzero friction case
2135 c
2136                ktm = 0.5d0*d_time/fricgam(i)
2137                psig = dsqrt(ktm*pterm) / fricgam(i)
2138                vsig = dsqrt(ktm*vterm)
2139                rhoc = dsqrt(1.0d0 - rho*rho)
2140                prand_vec(i) = psig
2141                vrand_vec1(i) = vsig * rho
2142                vrand_vec2(i) = vsig * rhoc
2143             end if
2144       end do
2145       if (lprn) then
2146       write (iout,*)
2147      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2148      &  " vrand_vec2"
2149       do i=1,dimen
2150         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2151      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2152       enddo
2153       endif
2154 c
2155 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2156 c
2157 #ifndef   LANG0
2158       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2159       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2160       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2161       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2162       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2163       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2164 #endif
2165 #ifdef MPI
2166       t_sdsetup=t_sdsetup+MPI_Wtime()
2167 #else
2168       t_sdsetup=t_sdsetup+tcpu()-tt0
2169 #endif
2170       return
2171       end
2172 c-------------------------------------------------------------
2173       subroutine eigtransf1(n,ndim,ab,d,c)
2174       implicit none
2175       integer n,ndim
2176       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2177       integer i,j,k
2178       do i=1,n
2179         do j=1,n
2180           c(i,j)=0.0d0
2181           do k=1,n
2182             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2183           enddo
2184         enddo
2185       enddo
2186       return
2187       end
2188 c-------------------------------------------------------------
2189       subroutine eigtransf(n,ndim,a,b,d,c)
2190       implicit none
2191       integer n,ndim
2192       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2193       integer i,j,k
2194       do i=1,n
2195         do j=1,n
2196           c(i,j)=0.0d0
2197           do k=1,n
2198             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2199           enddo
2200         enddo
2201       enddo
2202       return
2203       end
2204 c-------------------------------------------------------------
2205       subroutine sd_verlet1
2206 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2207       implicit real*8 (a-h,o-z)
2208       include 'DIMENSIONS'
2209       include 'COMMON.CONTROL'
2210       include 'COMMON.VAR'
2211       include 'COMMON.MD'
2212 #ifndef LANG0
2213       include 'COMMON.LANGEVIN'
2214 #else
2215       include 'COMMON.LANGEVIN.lang0'
2216 #endif
2217       include 'COMMON.CHAIN'
2218       include 'COMMON.DERIV'
2219       include 'COMMON.GEO'
2220       include 'COMMON.LOCAL'
2221       include 'COMMON.INTERACT'
2222       include 'COMMON.IOUNITS'
2223       include 'COMMON.NAMES'
2224       double precision stochforcvec(MAXRES6)
2225       common /stochcalc/ stochforcvec
2226       logical lprn /.false./
2227
2228 c      write (iout,*) "dc_old"
2229 c      do i=0,nres
2230 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2231 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2232 c      enddo
2233       do j=1,3
2234         dc_work(j)=dc_old(j,0)
2235         d_t_work(j)=d_t_old(j,0)
2236         d_a_work(j)=d_a_old(j,0)
2237       enddo
2238       ind=3
2239       do i=nnt,nct-1
2240         do j=1,3
2241           dc_work(ind+j)=dc_old(j,i)
2242           d_t_work(ind+j)=d_t_old(j,i)
2243           d_a_work(ind+j)=d_a_old(j,i)
2244         enddo
2245         ind=ind+3
2246       enddo
2247       do i=nnt,nct
2248         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2249           do j=1,3
2250             dc_work(ind+j)=dc_old(j,i+nres)
2251             d_t_work(ind+j)=d_t_old(j,i+nres)
2252             d_a_work(ind+j)=d_a_old(j,i+nres)
2253           enddo
2254           ind=ind+3
2255         endif
2256       enddo
2257 #ifndef LANG0
2258       if (lprn) then
2259       write (iout,*)
2260      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2261      &  " vrand_mat2"
2262       do i=1,dimen
2263         do j=1,dimen
2264           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2265      &      vfric_mat(i,j),afric_mat(i,j),
2266      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2267         enddo
2268       enddo
2269       endif
2270       do i=1,dimen
2271         ddt1=0.0d0
2272         ddt2=0.0d0
2273         do j=1,dimen
2274           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2275      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2276           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2277           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2278         enddo
2279         d_t_work_new(i)=ddt1+0.5d0*ddt2
2280         d_t_work(i)=ddt1+ddt2
2281       enddo
2282 #endif
2283       do j=1,3
2284         dc(j,0)=dc_work(j)
2285         d_t(j,0)=d_t_work(j)
2286       enddo
2287       ind=3
2288       do i=nnt,nct-1
2289         do j=1,3
2290           dc(j,i)=dc_work(ind+j)
2291           d_t(j,i)=d_t_work(ind+j)
2292         enddo
2293         ind=ind+3
2294       enddo
2295       do i=nnt,nct
2296         if (itype(i).ne.10) then
2297           inres=i+nres
2298           do j=1,3
2299             dc(j,inres)=dc_work(ind+j)
2300             d_t(j,inres)=d_t_work(ind+j)
2301           enddo
2302           ind=ind+3
2303         endif
2304       enddo
2305       return
2306       end
2307 c--------------------------------------------------------------------------
2308       subroutine sd_verlet2
2309 c  Calculating the adjusted velocities for accelerations
2310       implicit real*8 (a-h,o-z)
2311       include 'DIMENSIONS'
2312       include 'COMMON.CONTROL'
2313       include 'COMMON.VAR'
2314       include 'COMMON.MD'
2315 #ifndef LANG0
2316       include 'COMMON.LANGEVIN'
2317 #else
2318       include 'COMMON.LANGEVIN.lang0'
2319 #endif
2320       include 'COMMON.CHAIN'
2321       include 'COMMON.DERIV'
2322       include 'COMMON.GEO'
2323       include 'COMMON.LOCAL'
2324       include 'COMMON.INTERACT'
2325       include 'COMMON.IOUNITS'
2326       include 'COMMON.NAMES'
2327       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2328       common /stochcalc/ stochforcvec
2329 c
2330 c Compute the stochastic forces which contribute to velocity change
2331 c
2332       call stochastic_force(stochforcvecV)
2333
2334 #ifndef LANG0
2335       do i=1,dimen
2336         ddt1=0.0d0
2337         ddt2=0.0d0
2338         do j=1,dimen
2339           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2340           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2341      &     vrand_mat2(i,j)*stochforcvecV(j)
2342         enddo
2343         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2344       enddo
2345 #endif
2346       do j=1,3
2347         d_t(j,0)=d_t_work(j)
2348       enddo
2349       ind=3
2350       do i=nnt,nct-1
2351         do j=1,3
2352           d_t(j,i)=d_t_work(ind+j)
2353         enddo
2354         ind=ind+3
2355       enddo
2356       do i=nnt,nct
2357         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2358           inres=i+nres
2359           do j=1,3
2360             d_t(j,inres)=d_t_work(ind+j)
2361           enddo
2362           ind=ind+3
2363         endif
2364       enddo
2365       return
2366       end
2367 c-----------------------------------------------------------
2368       subroutine sd_verlet_ciccotti_setup
2369 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2370 c version
2371       implicit real*8 (a-h,o-z)
2372       include 'DIMENSIONS'
2373 #ifdef MPI
2374       include 'mpif.h'
2375 #endif
2376       include 'COMMON.CONTROL'
2377       include 'COMMON.VAR'
2378       include 'COMMON.MD'
2379 #ifndef LANG0
2380       include 'COMMON.LANGEVIN'
2381 #else
2382       include 'COMMON.LANGEVIN.lang0'
2383 #endif
2384       include 'COMMON.CHAIN'
2385       include 'COMMON.DERIV'
2386       include 'COMMON.GEO'
2387       include 'COMMON.LOCAL'
2388       include 'COMMON.INTERACT'
2389       include 'COMMON.IOUNITS'
2390       include 'COMMON.NAMES'
2391       include 'COMMON.TIME1'
2392       double precision emgdt(MAXRES6),
2393      & pterm,vterm,rho,rhoc,vsig,
2394      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2395      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2396      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2397       logical lprn /.false./
2398       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2399       double precision ktm
2400 #ifdef MPI
2401       tt0 = MPI_Wtime()
2402 #else
2403       tt0 = tcpu()
2404 #endif
2405 c
2406 c AL 8/17/04 Code adapted from tinker
2407 c
2408 c Get the frictional and random terms for stochastic dynamics in the
2409 c eigenspace of mass-scaled UNRES friction matrix
2410 c
2411       do i = 1, dimen
2412             write (iout,*) "i",i," fricgam",fricgam(i)
2413             gdt = fricgam(i) * d_time
2414 c
2415 c Stochastic dynamics reduces to simple MD for zero friction
2416 c
2417             if (gdt .le. zero) then
2418                pfric_vec(i) = 1.0d0
2419                vfric_vec(i) = d_time
2420                afric_vec(i) = 0.5d0*d_time*d_time
2421                prand_vec(i) = afric_vec(i)
2422                vrand_vec2(i) = vfric_vec(i)
2423 c
2424 c Analytical expressions when friction coefficient is large
2425 c
2426             else
2427                egdt = dexp(-gdt)
2428                pfric_vec(i) = egdt
2429                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2430                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2431                prand_vec(i) = afric_vec(i)
2432                vrand_vec2(i) = vfric_vec(i)
2433 c
2434 c Compute the scaling factors of random terms for the nonzero friction case
2435 c
2436 c               ktm = 0.5d0*d_time/fricgam(i)
2437 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2438 c               vsig = dsqrt(ktm*vterm)
2439 c               prand_vec(i) = psig*afric_vec(i)
2440 c               vrand_vec2(i) = vsig*vfric_vec(i)
2441             end if
2442       end do
2443       if (lprn) then
2444       write (iout,*)
2445      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2446      &  " vrand_vec2"
2447       do i=1,dimen
2448         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2449      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2450       enddo
2451       endif
2452 c
2453 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2454 c
2455       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2456       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2457       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2458       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2459       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2460 #ifdef MPI
2461       t_sdsetup=t_sdsetup+MPI_Wtime()
2462 #else
2463       t_sdsetup=t_sdsetup+tcpu()-tt0
2464 #endif
2465       return
2466       end
2467 c-------------------------------------------------------------
2468       subroutine sd_verlet1_ciccotti
2469 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2470       implicit real*8 (a-h,o-z)
2471       include 'DIMENSIONS'
2472 #ifdef MPI
2473       include 'mpif.h'
2474 #endif
2475       include 'COMMON.CONTROL'
2476       include 'COMMON.VAR'
2477       include 'COMMON.MD'
2478 #ifndef LANG0
2479       include 'COMMON.LANGEVIN'
2480 #else
2481       include 'COMMON.LANGEVIN.lang0'
2482 #endif
2483       include 'COMMON.CHAIN'
2484       include 'COMMON.DERIV'
2485       include 'COMMON.GEO'
2486       include 'COMMON.LOCAL'
2487       include 'COMMON.INTERACT'
2488       include 'COMMON.IOUNITS'
2489       include 'COMMON.NAMES'
2490       double precision stochforcvec(MAXRES6)
2491       common /stochcalc/ stochforcvec
2492       logical lprn /.false./
2493
2494 c      write (iout,*) "dc_old"
2495 c      do i=0,nres
2496 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2497 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2498 c      enddo
2499       do j=1,3
2500         dc_work(j)=dc_old(j,0)
2501         d_t_work(j)=d_t_old(j,0)
2502         d_a_work(j)=d_a_old(j,0)
2503       enddo
2504       ind=3
2505       do i=nnt,nct-1
2506         do j=1,3
2507           dc_work(ind+j)=dc_old(j,i)
2508           d_t_work(ind+j)=d_t_old(j,i)
2509           d_a_work(ind+j)=d_a_old(j,i)
2510         enddo
2511         ind=ind+3
2512       enddo
2513       do i=nnt,nct
2514         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2515           do j=1,3
2516             dc_work(ind+j)=dc_old(j,i+nres)
2517             d_t_work(ind+j)=d_t_old(j,i+nres)
2518             d_a_work(ind+j)=d_a_old(j,i+nres)
2519           enddo
2520           ind=ind+3
2521         endif
2522       enddo
2523
2524 #ifndef LANG0
2525       if (lprn) then
2526       write (iout,*)
2527      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2528      &  " vrand_mat2"
2529       do i=1,dimen
2530         do j=1,dimen
2531           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2532      &      vfric_mat(i,j),afric_mat(i,j),
2533      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2534         enddo
2535       enddo
2536       endif
2537       do i=1,dimen
2538         ddt1=0.0d0
2539         ddt2=0.0d0
2540         do j=1,dimen
2541           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2542      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2543           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2544           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2545         enddo
2546         d_t_work_new(i)=ddt1+0.5d0*ddt2
2547         d_t_work(i)=ddt1+ddt2
2548       enddo
2549 #endif
2550       do j=1,3
2551         dc(j,0)=dc_work(j)
2552         d_t(j,0)=d_t_work(j)
2553       enddo
2554       ind=3
2555       do i=nnt,nct-1
2556         do j=1,3
2557           dc(j,i)=dc_work(ind+j)
2558           d_t(j,i)=d_t_work(ind+j)
2559         enddo
2560         ind=ind+3
2561       enddo
2562       do i=nnt,nct
2563         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2564           inres=i+nres
2565           do j=1,3
2566             dc(j,inres)=dc_work(ind+j)
2567             d_t(j,inres)=d_t_work(ind+j)
2568           enddo
2569           ind=ind+3
2570         endif
2571       enddo
2572       return
2573       end
2574 c--------------------------------------------------------------------------
2575       subroutine sd_verlet2_ciccotti
2576 c  Calculating the adjusted velocities for accelerations
2577       implicit real*8 (a-h,o-z)
2578       include 'DIMENSIONS'
2579       include 'COMMON.CONTROL'
2580       include 'COMMON.VAR'
2581       include 'COMMON.MD'
2582 #ifndef LANG0
2583       include 'COMMON.LANGEVIN'
2584 #else
2585       include 'COMMON.LANGEVIN.lang0'
2586 #endif
2587       include 'COMMON.CHAIN'
2588       include 'COMMON.DERIV'
2589       include 'COMMON.GEO'
2590       include 'COMMON.LOCAL'
2591       include 'COMMON.INTERACT'
2592       include 'COMMON.IOUNITS'
2593       include 'COMMON.NAMES'
2594       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2595       common /stochcalc/ stochforcvec
2596 c
2597 c Compute the stochastic forces which contribute to velocity change
2598 c
2599       call stochastic_force(stochforcvecV)
2600 #ifndef LANG0
2601       do i=1,dimen
2602         ddt1=0.0d0
2603         ddt2=0.0d0
2604         do j=1,dimen
2605
2606           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2607 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2608           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2609         enddo
2610         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2611       enddo
2612 #endif
2613       do j=1,3
2614         d_t(j,0)=d_t_work(j)
2615       enddo
2616       ind=3
2617       do i=nnt,nct-1
2618         do j=1,3
2619           d_t(j,i)=d_t_work(ind+j)
2620         enddo
2621         ind=ind+3
2622       enddo
2623       do i=nnt,nct
2624         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2625           inres=i+nres
2626           do j=1,3
2627             d_t(j,inres)=d_t_work(ind+j)
2628           enddo
2629           ind=ind+3
2630         endif
2631       enddo
2632       return
2633       end
2634 #endif