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