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