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