source/unres/src-HCD-5D/MD_A-MTS.F

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