source/unres/src_MD-M/MD_A-MTS.F

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