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