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