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

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