source/unres/src_MD-M/MD.F

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