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