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          if (start_from_model) then
1613           i_model=iran_num(1,constr_homology)
1614           write (iout,*) 'starting from model ',i_model
1615           do i=1,2*nres
1616            do j=1,3
1617             c(j,i)=chomo(j,i,i_model)
1618            enddo
1619           enddo
1620
1621           call int_from_cart(.true.,.false.)
1622           call sc_loc_geom(.false.)
1623           do i=1,nres-1
1624             do j=1,3
1625              dc(j,i)=c(j,i+1)-c(j,i)
1626              dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1627             enddo
1628           enddo
1629           do i=2,nres-1
1630             do j=1,3
1631              dc(j,i+nres)=c(j,i+nres)-c(j,i)
1632              dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1633             enddo
1634           enddo
1635          endif
1636 c         call chainbuild_cart
1637          write (iout,*) "PREMINIM ",preminim
1638          if(iranconf.ne.0 .or. preminim) then
1639           if (overlapsc) then
1640            print *, 'Calling OVERLAP_SC'
1641            call overlap_sc(fail)
1642           endif
1643
1644           if (searchsc) then
1645            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1646            print *,'SC_move',nft_sc,etot
1647            if(me.eq.king.or..not.out1file)
1648      &      write(iout,*) 'SC_move',nft_sc,etot
1649           endif
1650
1651           if(dccart)then
1652            print *, 'Calling MINIM_DC'
1653            call minim_dc(etot,iretcode,nfun)
1654           else
1655            call geom_to_var(nvar,varia)
1656            print *,'Calling MINIMIZE.'
1657            call minimize(etot,varia,iretcode,nfun)
1658            call var_to_geom(nvar,varia)
1659           endif
1660           if(me.eq.king.or..not.out1file) then
1661              write(iout,*) "Minimized energy is",etot
1662              write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1663              call etotal(potEcomp)
1664              call enerprint(potEcomp)
1665           endif
1666          endif
1667       endif
1668       call chainbuild_cart
1669       call kinetic(EK)
1670       if (tbf) then
1671         call verlet_bath
1672       endif
1673       kinetic_T=2.0d0/(dimen3*Rb)*EK
1674       if(me.eq.king.or..not.out1file)then
1675        call cartprint
1676        call intout
1677       endif
1678 #ifdef MPI
1679       tt0=MPI_Wtime()
1680 #else
1681       tt0=tcpu()
1682 #endif
1683       call zerograd
1684       call etotal(potEcomp)
1685       if (large) call enerprint(potEcomp)
1686 #ifdef TIMING_ENE
1687 #ifdef MPI
1688       t_etotal=t_etotal+MPI_Wtime()-tt0
1689 #else
1690       t_etotal=t_etotal+tcpu()-tt0
1691 #endif
1692 #endif
1693       potE=potEcomp(0)
1694       call cartgrad
1695       call lagrangian
1696       call max_accel
1697       if (amax*d_time .gt. dvmax) then
1698         d_time=d_time*dvmax/amax
1699         if(me.eq.king.or..not.out1file) write (iout,*)
1700      &   "Time step reduced to",d_time,
1701      &   " because of too large initial acceleration."
1702       endif
1703 C      if(me.eq.king.or..not.out1file)then
1704 C       write(iout,*) "Potential energy and its components"
1705 C       call enerprint(potEcomp)
1706 c       write(iout,*) (potEcomp(i),i=0,n_ene)
1707 C      endif
1708       potE=potEcomp(0)-potEcomp(20)
1709       totE=EK+potE
1710       itime=0
1711       if (ntwe.ne.0) call statout(itime)
1712       if(me.eq.king.or..not.out1file)
1713      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
1714      &   " Kinetic energy",EK," Potential energy",potE,
1715      &   " Total energy",totE," Maximum acceleration ",
1716      &   amax
1717       if (large) then
1718         write (iout,*) "Initial coordinates"
1719         do i=1,nres
1720           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1721      &    (c(j,i+nres),j=1,3)
1722         enddo
1723         write (iout,*) "Initial dC"
1724         do i=0,nres
1725           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1726      &    (dc(j,i+nres),j=1,3)
1727         enddo
1728         write (iout,*) "Initial velocities"
1729         write (iout,"(13x,' backbone ',23x,' side chain')")
1730         do i=0,nres
1731           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1732      &    (d_t(j,i+nres),j=1,3)
1733         enddo
1734         write (iout,*) "Initial accelerations"
1735         do i=0,nres
1736 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1737           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
1738      &    (d_a(j,i+nres),j=1,3)
1739         enddo
1740       endif
1741       do i=0,2*nres
1742         do j=1,3
1743           dc_old(j,i)=dc(j,i)
1744           d_t_old(j,i)=d_t(j,i)
1745           d_a_old(j,i)=d_a(j,i)
1746         enddo
1747 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1748       enddo
1749       if (RESPA) then
1750 #ifdef MPI
1751         tt0 =MPI_Wtime()
1752 #else
1753         tt0 = tcpu()
1754 #endif
1755         call zerograd
1756         call etotal_short(energia_short)
1757         if (large) call enerprint(potEcomp)
1758 #ifdef TIMING_ENE
1759 #ifdef MPI
1760         t_eshort=t_eshort+MPI_Wtime()-tt0
1761 #else
1762         t_eshort=t_eshort+tcpu()-tt0
1763 #endif
1764 #endif
1765         call cartgrad
1766         call lagrangian
1767         if(.not.out1file .and. large) then
1768           write (iout,*) "energia_long",energia_long(0),
1769      &     " energia_short",energia_short(0),
1770      &     " total",energia_long(0)+energia_short(0)
1771           write (iout,*) "Initial fast-force accelerations"
1772           do i=0,nres
1773             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1774      &      (d_a(j,i+nres),j=1,3)
1775           enddo
1776         endif
1777 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1778         do i=0,2*nres
1779           do j=1,3
1780             d_a_short(j,i)=d_a(j,i)
1781           enddo
1782         enddo
1783 #ifdef MPI
1784         tt0=MPI_Wtime()
1785 #else
1786         tt0=tcpu()
1787 #endif
1788         call zerograd
1789         call etotal_long(energia_long)
1790         if (large) call enerprint(potEcomp)
1791 #ifdef TIMING_ENE
1792 #ifdef MPI
1793         t_elong=t_elong+MPI_Wtime()-tt0
1794 #else
1795         t_elong=t_elong+tcpu()-tt0
1796 #endif
1797 #endif
1798         call cartgrad
1799         call lagrangian
1800         if(.not.out1file .and. large) then
1801           write (iout,*) "energia_long",energia_long(0)
1802           write (iout,*) "Initial slow-force accelerations"
1803           do i=0,nres
1804             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1805      &      (d_a(j,i+nres),j=1,3)
1806           enddo
1807         endif
1808 #ifdef MPI
1809         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1810 #else
1811         t_enegrad=t_enegrad+tcpu()-tt0
1812 #endif
1813       endif
1814       return
1815       end
1816 c-----------------------------------------------------------
1817       subroutine random_vel
1818       implicit real*8 (a-h,o-z)
1819       include 'DIMENSIONS'
1820       include 'COMMON.CONTROL'
1821       include 'COMMON.VAR'
1822       include 'COMMON.MD'
1823 #ifndef LANG0
1824       include 'COMMON.LANGEVIN'
1825 #else
1826       include 'COMMON.LANGEVIN.lang0'
1827 #endif
1828       include 'COMMON.CHAIN'
1829       include 'COMMON.DERIV'
1830       include 'COMMON.GEO'
1831       include 'COMMON.LOCAL'
1832       include 'COMMON.INTERACT'
1833       include 'COMMON.IOUNITS'
1834       include 'COMMON.NAMES'
1835       include 'COMMON.TIME1'
1836       double precision xv,sigv,lowb,highb,vec_afm(3)
1837 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1838 c First generate velocities in the eigenspace of the G matrix
1839 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
1840 c      call flush(iout)
1841       xv=0.0d0
1842       ii=0
1843       do i=1,dimen
1844         do k=1,3
1845           ii=ii+1
1846           sigv=dsqrt((Rb*t_bath)/geigen(i))
1847           lowb=-5*sigv
1848           highb=5*sigv
1849           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
1850
1851 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
1852 c     &      " d_t_work_new",d_t_work_new(ii)
1853         enddo
1854       enddo
1855 C       if (SELFGUIDE.gt.0) then
1856 C       distance=0.0
1857 C       do j=1,3
1858 C       vec_afm(j)=c(j,afmend)-c(j,afmbeg)
1859 C       distance=distance+vec_afm(j)**2
1860 C       enddo
1861 C       distance=dsqrt(distance)
1862 C       do j=1,3
1863 C         d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
1864 C         d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
1865 C         write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
1866 C     &    d_t_work_new(j+(afmend-1)*3)
1867 C       enddo
1868
1869 C       endif
1870
1871 c diagnostics
1872 c      Ek1=0.0d0
1873 c      ii=0
1874 c      do i=1,dimen
1875 c        do k=1,3
1876 c          ii=ii+1
1877 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
1878 c        enddo
1879 c      enddo
1880 c      write (iout,*) "Ek from eigenvectors",Ek1
1881 c end diagnostics
1882 c Transform velocities to UNRES coordinate space
1883       do k=0,2
1884         do i=1,dimen
1885           ind=(i-1)*3+k+1
1886           d_t_work(ind)=0.0d0
1887           do j=1,dimen
1888             d_t_work(ind)=d_t_work(ind)
1889      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
1890           enddo
1891 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
1892 c          call flush(iout)
1893         enddo
1894       enddo
1895 c Transfer to the d_t vector
1896       do j=1,3
1897         d_t(j,0)=d_t_work(j)
1898       enddo
1899       ind=3
1900       do i=nnt,nct-1
1901         do j=1,3
1902           ind=ind+1
1903           d_t(j,i)=d_t_work(ind)
1904         enddo
1905       enddo
1906       do i=nnt,nct
1907         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1908           do j=1,3
1909             ind=ind+1
1910             d_t(j,i+nres)=d_t_work(ind)
1911           enddo
1912         endif
1913       enddo
1914 c      call kinetic(EK)
1915 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1916 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
1917 c      call flush(iout)
1918       return
1919       end
1920 #ifndef LANG0
1921 c-----------------------------------------------------------
1922       subroutine sd_verlet_p_setup
1923 c Sets up the parameters of stochastic Verlet algorithm
1924       implicit real*8 (a-h,o-z)
1925       include 'DIMENSIONS'
1926 #ifdef MPI
1927       include 'mpif.h'
1928 #endif
1929       include 'COMMON.CONTROL'
1930       include 'COMMON.VAR'
1931       include 'COMMON.MD'
1932 #ifndef LANG0
1933       include 'COMMON.LANGEVIN'
1934 #else
1935       include 'COMMON.LANGEVIN.lang0'
1936 #endif
1937       include 'COMMON.CHAIN'
1938       include 'COMMON.DERIV'
1939       include 'COMMON.GEO'
1940       include 'COMMON.LOCAL'
1941       include 'COMMON.INTERACT'
1942       include 'COMMON.IOUNITS'
1943       include 'COMMON.NAMES'
1944       include 'COMMON.TIME1'
1945       double precision emgdt(MAXRES6),
1946      & pterm,vterm,rho,rhoc,vsig,
1947      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1948      & afric_vec(MAXRES6),prand_vec(MAXRES6),
1949      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1950       logical lprn /.false./
1951       double precision zero /1.0d-8/, gdt_radius /0.05d0/
1952       double precision ktm
1953 #ifdef MPI
1954       tt0 = MPI_Wtime()
1955 #else
1956       tt0 = tcpu()
1957 #endif
1958 c
1959 c AL 8/17/04 Code adapted from tinker
1960 c
1961 c Get the frictional and random terms for stochastic dynamics in the
1962 c eigenspace of mass-scaled UNRES friction matrix
1963 c
1964       do i = 1, dimen
1965             gdt = fricgam(i) * d_time
1966 c
1967 c Stochastic dynamics reduces to simple MD for zero friction
1968 c
1969             if (gdt .le. zero) then
1970                pfric_vec(i) = 1.0d0
1971                vfric_vec(i) = d_time
1972                afric_vec(i) = 0.5d0 * d_time * d_time
1973                prand_vec(i) = 0.0d0
1974                vrand_vec1(i) = 0.0d0
1975                vrand_vec2(i) = 0.0d0
1976 c
1977 c Analytical expressions when friction coefficient is large
1978 c
1979             else
1980                if (gdt .ge. gdt_radius) then
1981                   egdt = dexp(-gdt)
1982                   pfric_vec(i) = egdt
1983                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
1984                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
1985                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
1986                   vterm = 1.0d0 - egdt**2
1987                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
1988 c
1989 c Use series expansions when friction coefficient is small
1990 c
1991                else
1992                   gdt2 = gdt * gdt
1993                   gdt3 = gdt * gdt2
1994                   gdt4 = gdt2 * gdt2
1995                   gdt5 = gdt2 * gdt3
1996                   gdt6 = gdt3 * gdt3
1997                   gdt7 = gdt3 * gdt4
1998                   gdt8 = gdt4 * gdt4
1999                   gdt9 = gdt4 * gdt5
2000                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2001      &                          - gdt5/120.0d0 + gdt6/720.0d0
2002      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2003      &                          - gdt9/362880.0d0) / fricgam(i)**2
2004                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2005                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2006                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2007      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2008      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2009      &                       + 127.0d0*gdt9/90720.0d0
2010                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2011      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2012      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2013      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2014                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2015      &                       - 17.0d0*gdt2/1280.0d0
2016      &                       + 17.0d0*gdt3/6144.0d0
2017      &                       + 40967.0d0*gdt4/34406400.0d0
2018      &                       - 57203.0d0*gdt5/275251200.0d0
2019      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2020                end if
2021 c
2022 c Compute the scaling factors of random terms for the nonzero friction case
2023 c
2024                ktm = 0.5d0*d_time/fricgam(i)
2025                psig = dsqrt(ktm*pterm) / fricgam(i)
2026                vsig = dsqrt(ktm*vterm)
2027                rhoc = dsqrt(1.0d0 - rho*rho)
2028                prand_vec(i) = psig
2029                vrand_vec1(i) = vsig * rho
2030                vrand_vec2(i) = vsig * rhoc
2031             end if
2032       end do
2033       if (lprn) then
2034       write (iout,*)
2035      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2036      &  " vrand_vec2"
2037       do i=1,dimen
2038         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2039      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2040       enddo
2041       endif
2042 c
2043 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2044 c
2045 #ifndef   LANG0
2046       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2047       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2048       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2049       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2050       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2051       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2052 #endif
2053 #ifdef MPI
2054       t_sdsetup=t_sdsetup+MPI_Wtime()
2055 #else
2056       t_sdsetup=t_sdsetup+tcpu()-tt0
2057 #endif
2058       return
2059       end
2060 c-------------------------------------------------------------
2061       subroutine eigtransf1(n,ndim,ab,d,c)
2062       implicit none
2063       integer n,ndim
2064       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2065       integer i,j,k
2066       do i=1,n
2067         do j=1,n
2068           c(i,j)=0.0d0
2069           do k=1,n
2070             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2071           enddo
2072         enddo
2073       enddo
2074       return
2075       end
2076 c-------------------------------------------------------------
2077       subroutine eigtransf(n,ndim,a,b,d,c)
2078       implicit none
2079       integer n,ndim
2080       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2081       integer i,j,k
2082       do i=1,n
2083         do j=1,n
2084           c(i,j)=0.0d0
2085           do k=1,n
2086             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2087           enddo
2088         enddo
2089       enddo
2090       return
2091       end
2092 c-------------------------------------------------------------
2093       subroutine sd_verlet1
2094 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2095       implicit real*8 (a-h,o-z)
2096       include 'DIMENSIONS'
2097       include 'COMMON.CONTROL'
2098       include 'COMMON.VAR'
2099       include 'COMMON.MD'
2100 #ifndef LANG0
2101       include 'COMMON.LANGEVIN'
2102 #else
2103       include 'COMMON.LANGEVIN.lang0'
2104 #endif
2105       include 'COMMON.CHAIN'
2106       include 'COMMON.DERIV'
2107       include 'COMMON.GEO'
2108       include 'COMMON.LOCAL'
2109       include 'COMMON.INTERACT'
2110       include 'COMMON.IOUNITS'
2111       include 'COMMON.NAMES'
2112       double precision stochforcvec(MAXRES6)
2113       common /stochcalc/ stochforcvec
2114       logical lprn /.false./
2115
2116 c      write (iout,*) "dc_old"
2117 c      do i=0,nres
2118 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2119 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2120 c      enddo
2121       do j=1,3
2122         dc_work(j)=dc_old(j,0)
2123         d_t_work(j)=d_t_old(j,0)
2124         d_a_work(j)=d_a_old(j,0)
2125       enddo
2126       ind=3
2127       do i=nnt,nct-1
2128         do j=1,3
2129           dc_work(ind+j)=dc_old(j,i)
2130           d_t_work(ind+j)=d_t_old(j,i)
2131           d_a_work(ind+j)=d_a_old(j,i)
2132         enddo
2133         ind=ind+3
2134       enddo
2135       do i=nnt,nct
2136         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2137           do j=1,3
2138             dc_work(ind+j)=dc_old(j,i+nres)
2139             d_t_work(ind+j)=d_t_old(j,i+nres)
2140             d_a_work(ind+j)=d_a_old(j,i+nres)
2141           enddo
2142           ind=ind+3
2143         endif
2144       enddo
2145 #ifndef LANG0
2146       if (lprn) then
2147       write (iout,*)
2148      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2149      &  " vrand_mat2"
2150       do i=1,dimen
2151         do j=1,dimen
2152           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2153      &      vfric_mat(i,j),afric_mat(i,j),
2154      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2155         enddo
2156       enddo
2157       endif
2158       do i=1,dimen
2159         ddt1=0.0d0
2160         ddt2=0.0d0
2161         do j=1,dimen
2162           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2163      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2164           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2165           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2166         enddo
2167         d_t_work_new(i)=ddt1+0.5d0*ddt2
2168         d_t_work(i)=ddt1+ddt2
2169       enddo
2170 #endif
2171       do j=1,3
2172         dc(j,0)=dc_work(j)
2173         d_t(j,0)=d_t_work(j)
2174       enddo
2175       ind=3
2176       do i=nnt,nct-1
2177         do j=1,3
2178           dc(j,i)=dc_work(ind+j)
2179           d_t(j,i)=d_t_work(ind+j)
2180         enddo
2181         ind=ind+3
2182       enddo
2183       do i=nnt,nct
2184         if (itype(i).ne.10) then
2185           inres=i+nres
2186           do j=1,3
2187             dc(j,inres)=dc_work(ind+j)
2188             d_t(j,inres)=d_t_work(ind+j)
2189           enddo
2190           ind=ind+3
2191         endif
2192       enddo
2193       return
2194       end
2195 c--------------------------------------------------------------------------
2196       subroutine sd_verlet2
2197 c  Calculating the adjusted velocities for accelerations
2198       implicit real*8 (a-h,o-z)
2199       include 'DIMENSIONS'
2200       include 'COMMON.CONTROL'
2201       include 'COMMON.VAR'
2202       include 'COMMON.MD'
2203 #ifndef LANG0
2204       include 'COMMON.LANGEVIN'
2205 #else
2206       include 'COMMON.LANGEVIN.lang0'
2207 #endif
2208       include 'COMMON.CHAIN'
2209       include 'COMMON.DERIV'
2210       include 'COMMON.GEO'
2211       include 'COMMON.LOCAL'
2212       include 'COMMON.INTERACT'
2213       include 'COMMON.IOUNITS'
2214       include 'COMMON.NAMES'
2215       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2216       common /stochcalc/ stochforcvec
2217 c
2218 c Compute the stochastic forces which contribute to velocity change
2219 c
2220       call stochastic_force(stochforcvecV)
2221
2222 #ifndef LANG0
2223       do i=1,dimen
2224         ddt1=0.0d0
2225         ddt2=0.0d0
2226         do j=1,dimen
2227           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2228           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2229      &     vrand_mat2(i,j)*stochforcvecV(j)
2230         enddo
2231         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2232       enddo
2233 #endif
2234       do j=1,3
2235         d_t(j,0)=d_t_work(j)
2236       enddo
2237       ind=3
2238       do i=nnt,nct-1
2239         do j=1,3
2240           d_t(j,i)=d_t_work(ind+j)
2241         enddo
2242         ind=ind+3
2243       enddo
2244       do i=nnt,nct
2245         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2246           inres=i+nres
2247           do j=1,3
2248             d_t(j,inres)=d_t_work(ind+j)
2249           enddo
2250           ind=ind+3
2251         endif
2252       enddo
2253       return
2254       end
2255 c-----------------------------------------------------------
2256       subroutine sd_verlet_ciccotti_setup
2257 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2258 c version
2259       implicit real*8 (a-h,o-z)
2260       include 'DIMENSIONS'
2261 #ifdef MPI
2262       include 'mpif.h'
2263 #endif
2264       include 'COMMON.CONTROL'
2265       include 'COMMON.VAR'
2266       include 'COMMON.MD'
2267 #ifndef LANG0
2268       include 'COMMON.LANGEVIN'
2269 #else
2270       include 'COMMON.LANGEVIN.lang0'
2271 #endif
2272       include 'COMMON.CHAIN'
2273       include 'COMMON.DERIV'
2274       include 'COMMON.GEO'
2275       include 'COMMON.LOCAL'
2276       include 'COMMON.INTERACT'
2277       include 'COMMON.IOUNITS'
2278       include 'COMMON.NAMES'
2279       include 'COMMON.TIME1'
2280       double precision emgdt(MAXRES6),
2281      & pterm,vterm,rho,rhoc,vsig,
2282      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2283      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2284      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2285       logical lprn /.false./
2286       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2287       double precision ktm
2288 #ifdef MPI
2289       tt0 = MPI_Wtime()
2290 #else
2291       tt0 = tcpu()
2292 #endif
2293 c
2294 c AL 8/17/04 Code adapted from tinker
2295 c
2296 c Get the frictional and random terms for stochastic dynamics in the
2297 c eigenspace of mass-scaled UNRES friction matrix
2298 c
2299       do i = 1, dimen
2300             write (iout,*) "i",i," fricgam",fricgam(i)
2301             gdt = fricgam(i) * d_time
2302 c
2303 c Stochastic dynamics reduces to simple MD for zero friction
2304 c
2305             if (gdt .le. zero) then
2306                pfric_vec(i) = 1.0d0
2307                vfric_vec(i) = d_time
2308                afric_vec(i) = 0.5d0*d_time*d_time
2309                prand_vec(i) = afric_vec(i)
2310                vrand_vec2(i) = vfric_vec(i)
2311 c
2312 c Analytical expressions when friction coefficient is large
2313 c
2314             else
2315                egdt = dexp(-gdt)
2316                pfric_vec(i) = egdt
2317                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2318                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2319                prand_vec(i) = afric_vec(i)
2320                vrand_vec2(i) = vfric_vec(i)
2321 c
2322 c Compute the scaling factors of random terms for the nonzero friction case
2323 c
2324 c               ktm = 0.5d0*d_time/fricgam(i)
2325 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2326 c               vsig = dsqrt(ktm*vterm)
2327 c               prand_vec(i) = psig*afric_vec(i)
2328 c               vrand_vec2(i) = vsig*vfric_vec(i)
2329             end if
2330       end do
2331       if (lprn) then
2332       write (iout,*)
2333      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2334      &  " vrand_vec2"
2335       do i=1,dimen
2336         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2337      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2338       enddo
2339       endif
2340 c
2341 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2342 c
2343       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2344       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2345       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2346       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2347       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2348 #ifdef MPI
2349       t_sdsetup=t_sdsetup+MPI_Wtime()
2350 #else
2351       t_sdsetup=t_sdsetup+tcpu()-tt0
2352 #endif
2353       return
2354       end
2355 c-------------------------------------------------------------
2356       subroutine sd_verlet1_ciccotti
2357 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2358       implicit real*8 (a-h,o-z)
2359       include 'DIMENSIONS'
2360 #ifdef MPI
2361       include 'mpif.h'
2362 #endif
2363       include 'COMMON.CONTROL'
2364       include 'COMMON.VAR'
2365       include 'COMMON.MD'
2366 #ifndef LANG0
2367       include 'COMMON.LANGEVIN'
2368 #else
2369       include 'COMMON.LANGEVIN.lang0'
2370 #endif
2371       include 'COMMON.CHAIN'
2372       include 'COMMON.DERIV'
2373       include 'COMMON.GEO'
2374       include 'COMMON.LOCAL'
2375       include 'COMMON.INTERACT'
2376       include 'COMMON.IOUNITS'
2377       include 'COMMON.NAMES'
2378       double precision stochforcvec(MAXRES6)
2379       common /stochcalc/ stochforcvec
2380       logical lprn /.false./
2381
2382 c      write (iout,*) "dc_old"
2383 c      do i=0,nres
2384 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2385 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2386 c      enddo
2387       do j=1,3
2388         dc_work(j)=dc_old(j,0)
2389         d_t_work(j)=d_t_old(j,0)
2390         d_a_work(j)=d_a_old(j,0)
2391       enddo
2392       ind=3
2393       do i=nnt,nct-1
2394         do j=1,3
2395           dc_work(ind+j)=dc_old(j,i)
2396           d_t_work(ind+j)=d_t_old(j,i)
2397           d_a_work(ind+j)=d_a_old(j,i)
2398         enddo
2399         ind=ind+3
2400       enddo
2401       do i=nnt,nct
2402         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2403           do j=1,3
2404             dc_work(ind+j)=dc_old(j,i+nres)
2405             d_t_work(ind+j)=d_t_old(j,i+nres)
2406             d_a_work(ind+j)=d_a_old(j,i+nres)
2407           enddo
2408           ind=ind+3
2409         endif
2410       enddo
2411
2412 #ifndef LANG0
2413       if (lprn) then
2414       write (iout,*)
2415      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2416      &  " vrand_mat2"
2417       do i=1,dimen
2418         do j=1,dimen
2419           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2420      &      vfric_mat(i,j),afric_mat(i,j),
2421      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2422         enddo
2423       enddo
2424       endif
2425       do i=1,dimen
2426         ddt1=0.0d0
2427         ddt2=0.0d0
2428         do j=1,dimen
2429           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2430      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2431           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2432           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2433         enddo
2434         d_t_work_new(i)=ddt1+0.5d0*ddt2
2435         d_t_work(i)=ddt1+ddt2
2436       enddo
2437 #endif
2438       do j=1,3
2439         dc(j,0)=dc_work(j)
2440         d_t(j,0)=d_t_work(j)
2441       enddo
2442       ind=3
2443       do i=nnt,nct-1
2444         do j=1,3
2445           dc(j,i)=dc_work(ind+j)
2446           d_t(j,i)=d_t_work(ind+j)
2447         enddo
2448         ind=ind+3
2449       enddo
2450       do i=nnt,nct
2451         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2452           inres=i+nres
2453           do j=1,3
2454             dc(j,inres)=dc_work(ind+j)
2455             d_t(j,inres)=d_t_work(ind+j)
2456           enddo
2457           ind=ind+3
2458         endif
2459       enddo
2460       return
2461       end
2462 c--------------------------------------------------------------------------
2463       subroutine sd_verlet2_ciccotti
2464 c  Calculating the adjusted velocities for accelerations
2465       implicit real*8 (a-h,o-z)
2466       include 'DIMENSIONS'
2467       include 'COMMON.CONTROL'
2468       include 'COMMON.VAR'
2469       include 'COMMON.MD'
2470 #ifndef LANG0
2471       include 'COMMON.LANGEVIN'
2472 #else
2473       include 'COMMON.LANGEVIN.lang0'
2474 #endif
2475       include 'COMMON.CHAIN'
2476       include 'COMMON.DERIV'
2477       include 'COMMON.GEO'
2478       include 'COMMON.LOCAL'
2479       include 'COMMON.INTERACT'
2480       include 'COMMON.IOUNITS'
2481       include 'COMMON.NAMES'
2482       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2483       common /stochcalc/ stochforcvec
2484 c
2485 c Compute the stochastic forces which contribute to velocity change
2486 c
2487       call stochastic_force(stochforcvecV)
2488 #ifndef LANG0
2489       do i=1,dimen
2490         ddt1=0.0d0
2491         ddt2=0.0d0
2492         do j=1,dimen
2493
2494           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2495 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2496           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2497         enddo
2498         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2499       enddo
2500 #endif
2501       do j=1,3
2502         d_t(j,0)=d_t_work(j)
2503       enddo
2504       ind=3
2505       do i=nnt,nct-1
2506         do j=1,3
2507           d_t(j,i)=d_t_work(ind+j)
2508         enddo
2509         ind=ind+3
2510       enddo
2511       do i=nnt,nct
2512         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2513           inres=i+nres
2514           do j=1,3
2515             d_t(j,inres)=d_t_work(ind+j)
2516           enddo
2517           ind=ind+3
2518         endif
2519       enddo
2520       return
2521       end
2522 #endif