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