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               if (me.eq.king.or..not.out1file) then
1969               write (iout,*) "Energies after removing overlaps"
1970               call etotal(energia(0))
1971               call enerprint(energia(0))
1972               endif
1973 #ifdef SEARCHSC
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           else
1991 ! Remove SC overlaps if requested
1992               if (overlapsc) then
1993                 write (iout,*) 'Calling OVERLAP_SC'
1994                 call overlap_sc(fail)
1995                 if (fail) then
1996                   write (iout,*)
1997      &            "Failed to remove overlap"
1998                 endif
1999               endif
2000               if (me.eq.king.or..not.out1file) then
2001               write (iout,*) "Energies after removing overlaps"
2002               call etotal(energia(0))
2003               call enerprint(energia(0))
2004               endif
2005           endif
2006 C 8/22/17 AL Minimize initial structure
2007           if (dccart) then
2008             if (me.eq.king.or..not.out1file) write(iout,*)
2009      &        'Minimizing initial PDB structure: Calling MINIM_DC'
2010             call minim_dc(etot,iretcode,nfun)
2011           else
2012             call geom_to_var(nvar,varia)
2013             if(me.eq.king.or..not.out1file) write (iout,*)
2014      &        'Minimizing initial PDB structure: Calling MINIMIZE.'
2015             call minimize(etot,varia,iretcode,nfun)
2016             call var_to_geom(nvar,varia)
2017 #ifdef LBFGS
2018             if (me.eq.king.or..not.out1file)
2019      &       write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2020             if(me.eq.king.or..not.out1file)
2021      &       write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2022 #else
2023             if (me.eq.king.or..not.out1file)
2024      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2025             if(me.eq.king.or..not.out1file)
2026      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2027 #endif
2028           endif
2029         endif
2030       endif ! .not. rest
2031       call chainbuild_cart
2032       call kinetic(EK)
2033       if (tbf) then
2034         call verlet_bath
2035       endif
2036       kinetic_T=2.0d0/(dimen3*Rb)*EK
2037       if(me.eq.king.or..not.out1file)then
2038        call cartprint
2039        call intout
2040       endif
2041 #ifdef MPI
2042       tt0=MPI_Wtime()
2043 #else
2044       tt0=tcpu()
2045 #endif
2046       call zerograd
2047       call etotal(potEcomp)
2048       if (large) call enerprint(potEcomp)
2049 #ifdef TIMING_ENE
2050 #ifdef MPI
2051       t_etotal=t_etotal+MPI_Wtime()-tt0
2052 #else
2053       t_etotal=t_etotal+tcpu()-tt0
2054 #endif
2055 #endif
2056       potE=potEcomp(0)
2057 c      write (iout,*) "PotE-homology",potE-potEcomp(27)
2058       call cartgrad
2059       call lagrangian
2060       call max_accel
2061       if (amax*d_time .gt. dvmax) then
2062         d_time=d_time*dvmax/amax
2063         if(me.eq.king.or..not.out1file) write (iout,*)
2064      &   "Time step reduced to",d_time,
2065      &   " because of too large initial acceleration."
2066       endif
2067       if(me.eq.king.or..not.out1file)then
2068        write(iout,*) "Potential energy and its components"
2069        call enerprint(potEcomp)
2070 c       write(iout,*) (potEcomp(i),i=0,n_ene)
2071       endif
2072       potE=potEcomp(0)-potEcomp(27)
2073 c      write (iout,*) "PotE-homology",potE
2074       totE=EK+potE
2075       itime=0
2076       itime_mat=itime
2077       if (ntwe.ne.0) call statout(itime)
2078       if(me.eq.king.or..not.out1file)
2079      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
2080      &   " Kinetic energy",EK," Potential energy",potE,
2081      &   " Total energy",totE," Maximum acceleration ",
2082      &   amax
2083       if (large) then
2084         write (iout,*) "Initial coordinates"
2085         do i=1,nres
2086           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
2087      &    (c(j,i+nres),j=1,3)
2088         enddo
2089         write (iout,*) "Initial dC"
2090         do i=0,nres
2091           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
2092      &    (dc(j,i+nres),j=1,3)
2093         enddo
2094         write (iout,*) "Initial velocities"
2095         write (iout,"(13x,' backbone ',23x,' side chain')")
2096         do i=0,nres
2097           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
2098      &    (d_t(j,i+nres),j=1,3)
2099         enddo
2100         write (iout,*) "Initial accelerations"
2101         do i=0,nres
2102 c          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2103           write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
2104      &    (d_a(j,i+nres),j=1,3)
2105         enddo
2106       endif
2107       do i=0,2*nres
2108         do j=1,3
2109           dc_old(j,i)=dc(j,i)
2110           d_t_old(j,i)=d_t(j,i)
2111           d_a_old(j,i)=d_a(j,i)
2112         enddo
2113 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2114       enddo
2115       if (RESPA) then
2116 #ifdef MPI
2117         tt0 =MPI_Wtime()
2118 #else
2119         tt0 = tcpu()
2120 #endif
2121         call zerograd
2122         call etotal_short(energia_short)
2123         if (large) call enerprint(potEcomp)
2124 #ifdef TIMING_ENE
2125 #ifdef MPI
2126         t_eshort=t_eshort+MPI_Wtime()-tt0
2127 #else
2128         t_eshort=t_eshort+tcpu()-tt0
2129 #endif
2130 #endif
2131         call cartgrad
2132         call lagrangian
2133         if(.not.out1file .and. large) then
2134           write (iout,*) "energia_long",energia_long(0),
2135      &     " energia_short",energia_short(0),
2136      &     " total",energia_long(0)+energia_short(0)
2137           write (iout,*) "Initial fast-force accelerations"
2138           do i=0,nres
2139             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2140      &      (d_a(j,i+nres),j=1,3)
2141           enddo
2142         endif
2143 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2144         do i=0,2*nres
2145           do j=1,3
2146             d_a_short(j,i)=d_a(j,i)
2147           enddo
2148         enddo
2149 #ifdef MPI
2150         tt0=MPI_Wtime()
2151 #else
2152         tt0=tcpu()
2153 #endif
2154         call zerograd
2155         call etotal_long(energia_long)
2156         if (large) call enerprint(potEcomp)
2157 #ifdef TIMING_ENE
2158 #ifdef MPI
2159         t_elong=t_elong+MPI_Wtime()-tt0
2160 #else
2161         t_elong=t_elong+tcpu()-tt0
2162 #endif
2163 #endif
2164         call cartgrad
2165         call lagrangian
2166         if(.not.out1file .and. large) then
2167           write (iout,*) "energia_long",energia_long(0)
2168           write (iout,*) "Initial slow-force accelerations"
2169           do i=0,nres
2170             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2171      &      (d_a(j,i+nres),j=1,3)
2172           enddo
2173         endif
2174 #ifdef MPI
2175         t_enegrad=t_enegrad+MPI_Wtime()-tt0
2176 #else
2177         t_enegrad=t_enegrad+tcpu()-tt0
2178 #endif
2179       endif
2180       return
2181       end
2182 c-----------------------------------------------------------
2183       subroutine random_vel
2184       implicit none
2185       include 'DIMENSIONS'
2186       include 'COMMON.CONTROL'
2187       include 'COMMON.VAR'
2188       include 'COMMON.MD'
2189 #ifdef FIVEDIAG
2190        include 'COMMON.LAGRANGE.5diag'
2191 #else
2192        include 'COMMON.LAGRANGE'
2193 #endif
2194 #ifndef LANG0
2195       include 'COMMON.LANGEVIN'
2196 #else
2197 #ifdef FIVEDIAG
2198       include 'COMMON.LANGEVIN.lang0.5diag'
2199 #else
2200       include 'COMMON.LANGEVIN.lang0'
2201 #endif
2202 #endif
2203       include 'COMMON.CHAIN'
2204       include 'COMMON.DERIV'
2205       include 'COMMON.GEO'
2206       include 'COMMON.LOCAL'
2207       include 'COMMON.INTERACT'
2208       include 'COMMON.IOUNITS'
2209       include 'COMMON.NAMES'
2210       include 'COMMON.TIME1'
2211       double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
2212       integer i,ii,j,k,l,ind
2213       double precision anorm_distr
2214       logical lprn /.false./
2215 #ifdef FIVEDIAG
2216       integer ichain,n,innt,inct,ibeg,ierr
2217       double precision work(8*maxres6)
2218       integer iwork(maxres6)
2219       double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
2220      & Gvec(maxres2_chain,maxres2_chain)
2221       common /przechowalnia/Ghalf,Geigen,Gvec
2222 #ifdef DEBUG
2223       double precision inertia(maxres2_chain,maxres2_chain)
2224 #endif
2225 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2226 c First generate velocities in the eigenspace of the G matrix
2227 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2228 c      call flush(iout)
2229 #ifdef DEBUG
2230       write (iout,*) "Random_vel, fivediag"
2231 #endif
2232       d_t=0.0d0
2233       Ek2=0.0d0
2234       EK=0.0d0
2235       Ek3=0.0d0
2236       do ichain=1,nchain
2237         ind=0
2238         ghalf=0.0d0
2239         n=dimen_chain(ichain)
2240         innt=iposd_chain(ichain)
2241         inct=innt+n-1
2242 #ifdef DEBUG
2243         write (iout,*) "Chain",ichain," n",n," start",innt
2244         do i=innt,inct
2245           if (i.lt.inct-1) then
2246             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
2247      &         DU2orig(i)
2248           else if (i.eq.inct-1) then
2249             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
2250           else
2251             write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
2252           endif
2253         enddo
2254 #endif
2255         ghalf(ind+1)=dmorig(innt)
2256         ghalf(ind+2)=du1orig(innt)
2257         ghalf(ind+3)=dmorig(innt+1)
2258         ind=ind+3
2259         do i=3,n
2260           ind=ind+i-3
2261 c          write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
2262 c     &       " indu1",innt+i-1," indm",innt+i
2263           ghalf(ind+1)=du2orig(innt-1+i-2)
2264           ghalf(ind+2)=du1orig(innt-1+i-1)
2265           ghalf(ind+3)=dmorig(innt-1+i)
2266 c          write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
2267 c     &       "DU1",innt-1+i-1,"DM ",innt-1+i
2268           ind=ind+3
2269         enddo
2270 #ifdef DEBUG
2271         ind=0
2272         do i=1,n
2273           do j=1,i
2274             ind=ind+1
2275             inertia(i,j)=ghalf(ind)
2276             inertia(j,i)=ghalf(ind)
2277           enddo
2278         enddo
2279 #endif
2280 #ifdef DEBUG
2281         write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
2282         write (iout,*) "Five-diagonal inertia matrix, lower triangle"
2283         call matoutr(n,ghalf)
2284 #endif
2285         call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
2286         if (large) then
2287           write (iout,'(//a,i3)')
2288      &    "Eigenvectors and eigenvalues of the G matrix chain",ichain
2289           call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
2290         endif
2291 #ifdef DIAGCHECK
2292 c check diagonalization
2293         do i=1,n
2294           do j=1,n
2295             aux=0.0d0
2296             do k=1,n
2297               do l=1,n
2298                 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
2299               enddo
2300             enddo
2301             if (i.eq.j) then
2302               write (iout,*) i,j,aux,geigen(i)
2303             else
2304               write (iout,*) i,j,aux
2305             endif
2306           enddo
2307         enddo
2308 #endif
2309         xv=0.0d0
2310         ii=0
2311         do i=1,n
2312           do k=1,3
2313             ii=ii+1
2314             sigv=dsqrt((Rb*t_bath)/geigen(i))
2315             lowb=-5*sigv
2316             highb=5*sigv
2317             d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2318             EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
2319 c            write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2320 c     &      " d_t_work_new",d_t_work_new(ii)
2321           enddo
2322         enddo
2323         do k=1,3
2324           do i=1,n
2325             ind=(i-1)*3+k
2326             d_t_work(ind)=0.0d0
2327             do j=1,n
2328               d_t_work(ind)=d_t_work(ind)
2329      &                +Gvec(i,j)*d_t_work_new((j-1)*3+k)
2330             enddo
2331 c            write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2332 c            call flush(iout)
2333           enddo
2334         enddo
2335 #ifdef DEBUG
2336         aux=0.0d0
2337         do k=1,3
2338           do i=1,n
2339             do j=1,n
2340             aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
2341             enddo
2342           enddo
2343         enddo
2344         Ek3=Ek3+aux/2
2345 #endif
2346 c Transfer to the d_t vector
2347         innt=chain_border(1,ichain)
2348         inct=chain_border(2,ichain)
2349         ind=0
2350 c        write (iout,*) "ichain",ichain," innt",innt," inct",inct
2351         do i=innt,inct
2352           do j=1,3
2353             ind=ind+1
2354             d_t(j,i)=d_t_work(ind)
2355           enddo
2356           if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2357             do j=1,3
2358               ind=ind+1
2359               d_t(j,i+nres)=d_t_work(ind)
2360             enddo
2361           endif
2362         enddo
2363       enddo
2364       if (large) then
2365         write (iout,*)
2366         write (iout,*) "Random velocities in the Calpha,SC space"
2367         do i=1,nres
2368           write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2369      &    restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2370         enddo
2371       endif
2372       call kinetic_CASC(Ek1)
2373 !
2374 ! Transform the velocities to virtual-bond space
2375 !
2376 #define WLOS
2377 #ifdef WLOS
2378       if (nnt.eq.1) then
2379         d_t(:,0)=d_t(:,1)
2380       endif
2381       do i=1,nres
2382         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
2383           do j=1,3
2384             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2385           enddo
2386         else
2387           do j=1,3
2388             d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2389             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2390           enddo
2391         end if
2392       enddo
2393       d_t(:,nres)=0.0d0
2394       d_t(:,nct)=0.0d0
2395       d_t(:,2*nres)=0.0d0
2396       if (nnt.gt.1) then
2397         d_t(:,0)=d_t(:,1)
2398         d_t(:,1)=0.0d0
2399       endif
2400 c      d_a(:,0)=d_a(:,1)
2401 c      d_a(:,1)=0.0d0
2402 c      write (iout,*) "Shifting accelerations"
2403       do ichain=2,nchain
2404 c        write (iout,*) "ichain",chain_border1(1,ichain)-1,
2405 c     &     chain_border1(1,ichain)
2406         d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
2407         d_t(:,chain_border1(1,ichain))=0.0d0
2408       enddo
2409 c      write (iout,*) "Adding accelerations"
2410       do ichain=2,nchain
2411 c        write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2412 c     &   chain_border(2,ichain-1)
2413         d_t(:,chain_border1(1,ichain)-1)=
2414      &  d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2415         d_t(:,chain_border(2,ichain-1))=0.0d0
2416       enddo
2417       do ichain=2,nchain
2418         write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
2419      &   chain_border(2,ichain-1)
2420         d_t(:,chain_border1(1,ichain)-1)=
2421      &  d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
2422         d_t(:,chain_border(2,ichain-1))=0.0d0
2423       enddo
2424 #else
2425       ibeg=0
2426 c      do j=1,3
2427 c        d_t(j,0)=d_t(j,nnt)
2428 c      enddo
2429       do ichain=1,nchain
2430       innt=chain_border(1,ichain)
2431       inct=chain_border(2,ichain)
2432 c      write (iout,*) "ichain",ichain," innt",innt," inct",inct
2433 c      write (iout,*) "ibeg",ibeg
2434       do j=1,3
2435         d_t(j,ibeg)=d_t(j,innt)
2436       enddo
2437       ibeg=inct+1
2438       do i=innt,inct
2439         if (iabs(itype(i).eq.10)) then
2440 c          write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
2441           do j=1,3
2442             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2443           enddo
2444         else
2445           do j=1,3
2446             d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2447             d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2448           enddo
2449         end if
2450       enddo
2451       enddo
2452 #endif
2453       if (large) then
2454         write (iout,*)
2455         write (iout,*)
2456      &    "Random velocities in the virtual-bond-vector space"
2457         write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
2458         do i=1,nres
2459           write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
2460      &     restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
2461         enddo
2462         write (iout,*)
2463         write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
2464      &   Ek
2465         write (iout,*)
2466      &   "Kinetic temperatures from inertia matrix eigenvalues",
2467      &   2*Ek/(3*dimen*Rb)
2468 #ifdef DEBUG
2469         write (iout,*) "Kinetic energy from inertia matrix",Ek3
2470         write (iout,*) "Kinetic temperatures from inertia",
2471      &   2*Ek3/(3*dimen*Rb)
2472 #endif
2473         write (iout,*) "Kinetic energy from velocities in CA-SC space",
2474      &   Ek1
2475         write (iout,*)
2476      &   "Kinetic temperatures from velovities in CA-SC space",
2477      &   2*Ek1/(3*dimen*Rb)
2478         call kinetic(Ek1)
2479         write (iout,*)
2480      &   "Kinetic energy from virtual-bond-vector velocities",Ek1
2481         write (iout,*)
2482      &   "Kinetic temperature from virtual-bond-vector velocities ",
2483      &   2*Ek1/(dimen3*Rb)
2484       endif
2485 #else
2486       xv=0.0d0
2487       ii=0
2488       do i=1,dimen
2489         do k=1,3
2490           ii=ii+1
2491           sigv=dsqrt((Rb*t_bath)/geigen(i))
2492           lowb=-5*sigv
2493           highb=5*sigv
2494           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2495
2496 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
2497 c     &      " d_t_work_new",d_t_work_new(ii)
2498         enddo
2499       enddo
2500 C       if (SELFGUIDE.gt.0) then
2501 C       distance=0.0
2502 C       do j=1,3
2503 C       vec_afm(j)=c(j,afmend)-c(j,afmbeg)
2504 C       distance=distance+vec_afm(j)**2
2505 C       enddo
2506 C       distance=dsqrt(distance)
2507 C       do j=1,3
2508 C         d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
2509 C         d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
2510 C         write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
2511 C     &    d_t_work_new(j+(afmend-1)*3)
2512 C       enddo
2513
2514 C       endif
2515
2516 c diagnostics
2517 c      Ek1=0.0d0
2518 c      ii=0
2519 c      do i=1,dimen
2520 c        do k=1,3
2521 c          ii=ii+1
2522 c          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2523 c        enddo
2524 c      enddo
2525 c      write (iout,*) "Ek from eigenvectors",Ek1
2526 c end diagnostics
2527 c Transform velocities to UNRES coordinate space
2528       do k=0,2
2529         do i=1,dimen
2530           ind=(i-1)*3+k+1
2531           d_t_work(ind)=0.0d0
2532           do j=1,dimen
2533             d_t_work(ind)=d_t_work(ind)
2534      &                      +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2535           enddo
2536 c          write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2537 c          call flush(iout)
2538         enddo
2539       enddo
2540 c Transfer to the d_t vector
2541       do j=1,3
2542         d_t(j,0)=d_t_work(j)
2543       enddo
2544       ind=3
2545       do i=nnt,nct-1
2546         do j=1,3
2547           ind=ind+1
2548           d_t(j,i)=d_t_work(ind)
2549         enddo
2550       enddo
2551       do i=nnt,nct
2552         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2553           do j=1,3
2554             ind=ind+1
2555             d_t(j,i+nres)=d_t_work(ind)
2556           enddo
2557         endif
2558       enddo
2559 c      call kinetic(EK)
2560 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
2561 c     &  2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2562 c      call flush(iout)
2563 #endif
2564       return
2565       end
2566 #ifndef LANG0
2567 c-----------------------------------------------------------
2568       subroutine sd_verlet_p_setup
2569 c Sets up the parameters of stochastic Verlet algorithm
2570       implicit none
2571       include 'DIMENSIONS'
2572 #ifdef MPI
2573       include 'mpif.h'
2574 #endif
2575       include 'COMMON.CONTROL'
2576       include 'COMMON.VAR'
2577       include 'COMMON.MD'
2578 #ifndef LANG0
2579       include 'COMMON.LANGEVIN'
2580 #else
2581 #ifdef FIVEDIAG
2582       include 'COMMON.LANGEVIN.lang0.5diag'
2583 #else
2584       include 'COMMON.LANGEVIN.lang0'
2585 #endif
2586 #endif
2587       include 'COMMON.CHAIN'
2588       include 'COMMON.DERIV'
2589       include 'COMMON.GEO'
2590       include 'COMMON.LOCAL'
2591       include 'COMMON.INTERACT'
2592       include 'COMMON.IOUNITS'
2593       include 'COMMON.NAMES'
2594       include 'COMMON.TIME1'
2595       double precision emgdt(MAXRES6),
2596      & pterm,vterm,rho,rhoc,vsig,
2597      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2598      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2599      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2600       logical lprn /.false./
2601       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2602       double precision ktm
2603 #ifdef MPI
2604       tt0 = MPI_Wtime()
2605 #else
2606       tt0 = tcpu()
2607 #endif
2608 c
2609 c AL 8/17/04 Code adapted from tinker
2610 c
2611 c Get the frictional and random terms for stochastic dynamics in the
2612 c eigenspace of mass-scaled UNRES friction matrix
2613 c
2614       do i = 1, dimen
2615             gdt = fricgam(i) * d_time
2616 c
2617 c Stochastic dynamics reduces to simple MD for zero friction
2618 c
2619             if (gdt .le. zero) then
2620                pfric_vec(i) = 1.0d0
2621                vfric_vec(i) = d_time
2622                afric_vec(i) = 0.5d0 * d_time * d_time
2623                prand_vec(i) = 0.0d0
2624                vrand_vec1(i) = 0.0d0
2625                vrand_vec2(i) = 0.0d0
2626 c
2627 c Analytical expressions when friction coefficient is large
2628 c
2629             else
2630                if (gdt .ge. gdt_radius) then
2631                   egdt = dexp(-gdt)
2632                   pfric_vec(i) = egdt
2633                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2634                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2635                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2636                   vterm = 1.0d0 - egdt**2
2637                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2638 c
2639 c Use series expansions when friction coefficient is small
2640 c
2641                else
2642                   gdt2 = gdt * gdt
2643                   gdt3 = gdt * gdt2
2644                   gdt4 = gdt2 * gdt2
2645                   gdt5 = gdt2 * gdt3
2646                   gdt6 = gdt3 * gdt3
2647                   gdt7 = gdt3 * gdt4
2648                   gdt8 = gdt4 * gdt4
2649                   gdt9 = gdt4 * gdt5
2650                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2651      &                          - gdt5/120.0d0 + gdt6/720.0d0
2652      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2653      &                          - gdt9/362880.0d0) / fricgam(i)**2
2654                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2655                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2656                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2657      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2658      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2659      &                       + 127.0d0*gdt9/90720.0d0
2660                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2661      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2662      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2663      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2664                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2665      &                       - 17.0d0*gdt2/1280.0d0
2666      &                       + 17.0d0*gdt3/6144.0d0
2667      &                       + 40967.0d0*gdt4/34406400.0d0
2668      &                       - 57203.0d0*gdt5/275251200.0d0
2669      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2670                end if
2671 c
2672 c Compute the scaling factors of random terms for the nonzero friction case
2673 c
2674                ktm = 0.5d0*d_time/fricgam(i)
2675                psig = dsqrt(ktm*pterm) / fricgam(i)
2676                vsig = dsqrt(ktm*vterm)
2677                rhoc = dsqrt(1.0d0 - rho*rho)
2678                prand_vec(i) = psig
2679                vrand_vec1(i) = vsig * rho
2680                vrand_vec2(i) = vsig * rhoc
2681             end if
2682       end do
2683       if (lprn) then
2684       write (iout,*)
2685      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2686      &  " vrand_vec2"
2687       do i=1,dimen
2688         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2689      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2690       enddo
2691       endif
2692 c
2693 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2694 c
2695       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2696       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2697       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2698       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2699       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2700       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2701 #ifdef MPI
2702       t_sdsetup=t_sdsetup+MPI_Wtime()
2703 #else
2704       t_sdsetup=t_sdsetup+tcpu()-tt0
2705 #endif
2706       return
2707       end
2708 c-------------------------------------------------------------
2709       subroutine eigtransf1(n,ndim,ab,d,c)
2710       implicit none
2711       integer n,ndim
2712       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2713       integer i,j,k
2714       do i=1,n
2715         do j=1,n
2716           c(i,j)=0.0d0
2717           do k=1,n
2718             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2719           enddo
2720         enddo
2721       enddo
2722       return
2723       end
2724 c-------------------------------------------------------------
2725       subroutine eigtransf(n,ndim,a,b,d,c)
2726       implicit none
2727       integer n,ndim
2728       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2729       integer i,j,k
2730       do i=1,n
2731         do j=1,n
2732           c(i,j)=0.0d0
2733           do k=1,n
2734             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2735           enddo
2736         enddo
2737       enddo
2738       return
2739       end
2740 c-------------------------------------------------------------
2741       subroutine sd_verlet1
2742 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2743       implicit none
2744       include 'DIMENSIONS'
2745       include 'COMMON.CONTROL'
2746       include 'COMMON.VAR'
2747       include 'COMMON.MD'
2748 #ifdef FIVEDIAG
2749        include 'COMMON.LAGRANGE.5diag'
2750 #else
2751        include 'COMMON.LAGRANGE'
2752 #endif
2753 #ifndef LANG0
2754       include 'COMMON.LANGEVIN'
2755 #else
2756 #ifdef FIVEDIAG
2757       include 'COMMON.LANGEVIN.lang0.5diag'
2758 #else
2759       include 'COMMON.LANGEVIN.lang0'
2760 #endif
2761 #endif
2762       include 'COMMON.CHAIN'
2763       include 'COMMON.DERIV'
2764       include 'COMMON.GEO'
2765       include 'COMMON.LOCAL'
2766       include 'COMMON.INTERACT'
2767       include 'COMMON.IOUNITS'
2768       include 'COMMON.NAMES'
2769       double precision stochforcvec(MAXRES6)
2770       common /stochcalc/ stochforcvec
2771       logical lprn /.false./
2772       integer i,j,ind,inres
2773
2774 c      write (iout,*) "dc_old"
2775 c      do i=0,nres
2776 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2777 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2778 c      enddo
2779       do j=1,3
2780         dc_work(j)=dc_old(j,0)
2781         d_t_work(j)=d_t_old(j,0)
2782         d_a_work(j)=d_a_old(j,0)
2783       enddo
2784       ind=3
2785       do i=nnt,nct-1
2786         do j=1,3
2787           dc_work(ind+j)=dc_old(j,i)
2788           d_t_work(ind+j)=d_t_old(j,i)
2789           d_a_work(ind+j)=d_a_old(j,i)
2790         enddo
2791         ind=ind+3
2792       enddo
2793       do i=nnt,nct
2794         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2795           do j=1,3
2796             dc_work(ind+j)=dc_old(j,i+nres)
2797             d_t_work(ind+j)=d_t_old(j,i+nres)
2798             d_a_work(ind+j)=d_a_old(j,i+nres)
2799           enddo
2800           ind=ind+3
2801         endif
2802       enddo
2803 #ifndef LANG0
2804       if (lprn) then
2805       write (iout,*)
2806      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2807      &  " vrand_mat2"
2808       do i=1,dimen
2809         do j=1,dimen
2810           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2811      &      vfric_mat(i,j),afric_mat(i,j),
2812      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2813         enddo
2814       enddo
2815       endif
2816       do i=1,dimen
2817         ddt1=0.0d0
2818         ddt2=0.0d0
2819         do j=1,dimen
2820           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2821      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2822           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2823           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2824         enddo
2825         d_t_work_new(i)=ddt1+0.5d0*ddt2
2826         d_t_work(i)=ddt1+ddt2
2827       enddo
2828 #endif
2829       do j=1,3
2830         dc(j,0)=dc_work(j)
2831         d_t(j,0)=d_t_work(j)
2832       enddo
2833       ind=3
2834       do i=nnt,nct-1
2835         do j=1,3
2836           dc(j,i)=dc_work(ind+j)
2837           d_t(j,i)=d_t_work(ind+j)
2838         enddo
2839         ind=ind+3
2840       enddo
2841       do i=nnt,nct
2842         if (itype(i).ne.10) then
2843           inres=i+nres
2844           do j=1,3
2845             dc(j,inres)=dc_work(ind+j)
2846             d_t(j,inres)=d_t_work(ind+j)
2847           enddo
2848           ind=ind+3
2849         endif
2850       enddo
2851       return
2852       end
2853 c--------------------------------------------------------------------------
2854       subroutine sd_verlet2
2855 c  Calculating the adjusted velocities for accelerations
2856       implicit none
2857       include 'DIMENSIONS'
2858       include 'COMMON.CONTROL'
2859       include 'COMMON.VAR'
2860       include 'COMMON.MD'
2861 #ifdef FIVEDIAG
2862        include 'COMMON.LAGRANGE.5diag'
2863 #else
2864        include 'COMMON.LAGRANGE'
2865 #endif
2866 #ifndef LANG0
2867       include 'COMMON.LANGEVIN'
2868 #else
2869 #ifdef FIVEDIAG
2870       include 'COMMON.LANGEVIN.lang0.5diag'
2871 #else
2872       include 'COMMON.LANGEVIN.lang0'
2873 #endif
2874 #endif
2875       include 'COMMON.CHAIN'
2876       include 'COMMON.DERIV'
2877       include 'COMMON.GEO'
2878       include 'COMMON.LOCAL'
2879       include 'COMMON.INTERACT'
2880       include 'COMMON.IOUNITS'
2881       include 'COMMON.NAMES'
2882       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2883       common /stochcalc/ stochforcvec
2884       integer i,j,ind,inres
2885 c
2886 c Compute the stochastic forces which contribute to velocity change
2887 c
2888       call stochastic_force(stochforcvecV)
2889
2890 #ifndef LANG0
2891       do i=1,dimen
2892         ddt1=0.0d0
2893         ddt2=0.0d0
2894         do j=1,dimen
2895           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2896           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2897      &     vrand_mat2(i,j)*stochforcvecV(j)
2898         enddo
2899         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2900       enddo
2901 #endif
2902       do j=1,3
2903         d_t(j,0)=d_t_work(j)
2904       enddo
2905       ind=3
2906       do i=nnt,nct-1
2907         do j=1,3
2908           d_t(j,i)=d_t_work(ind+j)
2909         enddo
2910         ind=ind+3
2911       enddo
2912       do i=nnt,nct
2913         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2914           inres=i+nres
2915           do j=1,3
2916             d_t(j,inres)=d_t_work(ind+j)
2917           enddo
2918           ind=ind+3
2919         endif
2920       enddo
2921       return
2922       end
2923 c-----------------------------------------------------------
2924       subroutine sd_verlet_ciccotti_setup
2925 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2926 c version
2927       implicit none
2928       include 'DIMENSIONS'
2929 #ifdef MPI
2930       include 'mpif.h'
2931 #endif
2932       include 'COMMON.CONTROL'
2933       include 'COMMON.VAR'
2934       include 'COMMON.MD'
2935 #ifdef FIVEDIAG
2936       include 'COMMON.LAGRANGE.5diag'
2937 #else
2938       include 'COMMON.LAGRANGE'
2939 #endif
2940       include 'COMMON.LANGEVIN'
2941       include 'COMMON.CHAIN'
2942       include 'COMMON.DERIV'
2943       include 'COMMON.GEO'
2944       include 'COMMON.LOCAL'
2945       include 'COMMON.INTERACT'
2946       include 'COMMON.IOUNITS'
2947       include 'COMMON.NAMES'
2948       include 'COMMON.TIME1'
2949       double precision emgdt(MAXRES6),
2950      & pterm,vterm,rho,rhoc,vsig,
2951      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2952      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2953      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2954       logical lprn /.false./
2955       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2956       double precision ktm
2957       integer i
2958 #ifdef MPI
2959       tt0 = MPI_Wtime()
2960 #else
2961       tt0 = tcpu()
2962 #endif
2963 c
2964 c AL 8/17/04 Code adapted from tinker
2965 c
2966 c Get the frictional and random terms for stochastic dynamics in the
2967 c eigenspace of mass-scaled UNRES friction matrix
2968 c
2969       do i = 1, dimen
2970             write (iout,*) "i",i," fricgam",fricgam(i)
2971             gdt = fricgam(i) * d_time
2972 c
2973 c Stochastic dynamics reduces to simple MD for zero friction
2974 c
2975             if (gdt .le. zero) then
2976                pfric_vec(i) = 1.0d0
2977                vfric_vec(i) = d_time
2978                afric_vec(i) = 0.5d0*d_time*d_time
2979                prand_vec(i) = afric_vec(i)
2980                vrand_vec2(i) = vfric_vec(i)
2981 c
2982 c Analytical expressions when friction coefficient is large
2983 c
2984             else
2985                egdt = dexp(-gdt)
2986                pfric_vec(i) = egdt
2987                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2988                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2989                prand_vec(i) = afric_vec(i)
2990                vrand_vec2(i) = vfric_vec(i)
2991 c
2992 c Compute the scaling factors of random terms for the nonzero friction case
2993 c
2994 c               ktm = 0.5d0*d_time/fricgam(i)
2995 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2996 c               vsig = dsqrt(ktm*vterm)
2997 c               prand_vec(i) = psig*afric_vec(i)
2998 c               vrand_vec2(i) = vsig*vfric_vec(i)
2999             end if
3000       end do
3001       if (lprn) then
3002       write (iout,*)
3003      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
3004      &  " vrand_vec2"
3005       do i=1,dimen
3006         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
3007      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3008       enddo
3009       endif
3010 c
3011 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3012 c
3013       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3014       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3015       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3016       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3017       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3018 #ifdef MPI
3019       t_sdsetup=t_sdsetup+MPI_Wtime()
3020 #else
3021       t_sdsetup=t_sdsetup+tcpu()-tt0
3022 #endif
3023       return
3024       end
3025 c-------------------------------------------------------------
3026       subroutine sd_verlet1_ciccotti
3027 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
3028       implicit none
3029       include 'DIMENSIONS'
3030 #ifdef MPI
3031       include 'mpif.h'
3032 #endif
3033       include 'COMMON.CONTROL'
3034       include 'COMMON.VAR'
3035       include 'COMMON.MD'
3036 #ifdef FIVEDIAG
3037       include 'COMMON.LAGRANGE.5diag'
3038 #else
3039       include 'COMMON.LAGRANGE'
3040 #endif
3041       include 'COMMON.LANGEVIN'
3042       include 'COMMON.CHAIN'
3043       include 'COMMON.DERIV'
3044       include 'COMMON.GEO'
3045       include 'COMMON.LOCAL'
3046       include 'COMMON.INTERACT'
3047       include 'COMMON.IOUNITS'
3048       include 'COMMON.NAMES'
3049       double precision stochforcvec(MAXRES6)
3050       common /stochcalc/ stochforcvec
3051       logical lprn /.false./
3052       integer i,j
3053
3054 c      write (iout,*) "dc_old"
3055 c      do i=0,nres
3056 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
3057 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3058 c      enddo
3059       do j=1,3
3060         dc_work(j)=dc_old(j,0)
3061         d_t_work(j)=d_t_old(j,0)
3062         d_a_work(j)=d_a_old(j,0)
3063       enddo
3064       ind=3
3065       do i=nnt,nct-1
3066         do j=1,3
3067           dc_work(ind+j)=dc_old(j,i)
3068           d_t_work(ind+j)=d_t_old(j,i)
3069           d_a_work(ind+j)=d_a_old(j,i)
3070         enddo
3071         ind=ind+3
3072       enddo
3073       do i=nnt,nct
3074         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3075           do j=1,3
3076             dc_work(ind+j)=dc_old(j,i+nres)
3077             d_t_work(ind+j)=d_t_old(j,i+nres)
3078             d_a_work(ind+j)=d_a_old(j,i+nres)
3079           enddo
3080           ind=ind+3
3081         endif
3082       enddo
3083
3084 #ifndef LANG0
3085       if (lprn) then
3086       write (iout,*)
3087      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
3088      &  " vrand_mat2"
3089       do i=1,dimen
3090         do j=1,dimen
3091           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
3092      &      vfric_mat(i,j),afric_mat(i,j),
3093      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3094         enddo
3095       enddo
3096       endif
3097       do i=1,dimen
3098         ddt1=0.0d0
3099         ddt2=0.0d0
3100         do j=1,dimen
3101           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
3102      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3103           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3104           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3105         enddo
3106         d_t_work_new(i)=ddt1+0.5d0*ddt2
3107         d_t_work(i)=ddt1+ddt2
3108       enddo
3109 #endif
3110       do j=1,3
3111         dc(j,0)=dc_work(j)
3112         d_t(j,0)=d_t_work(j)
3113       enddo
3114       ind=3
3115       do i=nnt,nct-1
3116         do j=1,3
3117           dc(j,i)=dc_work(ind+j)
3118           d_t(j,i)=d_t_work(ind+j)
3119         enddo
3120         ind=ind+3
3121       enddo
3122       do i=nnt,nct
3123         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3124           inres=i+nres
3125           do j=1,3
3126             dc(j,inres)=dc_work(ind+j)
3127             d_t(j,inres)=d_t_work(ind+j)
3128           enddo
3129           ind=ind+3
3130         endif
3131       enddo
3132       return
3133       end
3134 c--------------------------------------------------------------------------
3135       subroutine sd_verlet2_ciccotti
3136 c  Calculating the adjusted velocities for accelerations
3137       implicit none
3138       include 'DIMENSIONS'
3139       include 'COMMON.CONTROL'
3140       include 'COMMON.VAR'
3141       include 'COMMON.MD'
3142 #ifdef FIVEDIAG
3143       include 'COMMON.LAGRANGE.5diag'
3144 #else
3145       include 'COMMON.LAGRANGE'
3146 #endif
3147       include 'COMMON.LANGEVIN'
3148       include 'COMMON.CHAIN'
3149       include 'COMMON.DERIV'
3150       include 'COMMON.GEO'
3151       include 'COMMON.LOCAL'
3152       include 'COMMON.INTERACT'
3153       include 'COMMON.IOUNITS'
3154       include 'COMMON.NAMES'
3155       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
3156       common /stochcalc/ stochforcvec
3157       integer i,j
3158 c
3159 c Compute the stochastic forces which contribute to velocity change
3160 c
3161       call stochastic_force(stochforcvecV)
3162 #ifndef LANG0
3163       do i=1,dimen
3164         ddt1=0.0d0
3165         ddt2=0.0d0
3166         do j=1,dimen
3167
3168           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3169 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3170           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3171         enddo
3172         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3173       enddo
3174 #endif
3175       do j=1,3
3176         d_t(j,0)=d_t_work(j)
3177       enddo
3178       ind=3
3179       do i=nnt,nct-1
3180         do j=1,3
3181           d_t(j,i)=d_t_work(ind+j)
3182         enddo
3183         ind=ind+3
3184       enddo
3185       do i=nnt,nct
3186         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3187           inres=i+nres
3188           do j=1,3
3189             d_t(j,inres)=d_t_work(ind+j)
3190           enddo
3191           ind=ind+3
3192         endif
3193       enddo
3194       return
3195       end
3196 #endif