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           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(iabs(itype(i))))
1635         enddo
1636       endif
1637 c Open the pdb file for snapshotshots
1638 #ifdef MPI
1639       if(mdpdb) then
1640         if (ilen(tmpdir).gt.0)
1641      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1642      &      liczba(:ilen(liczba))//".pdb")
1643         open(ipdb,
1644      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1645      &  //".pdb")
1646       else
1647 #ifdef NOXDR
1648         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1649      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1650      &      liczba(:ilen(liczba))//".x")
1651         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1652      &  //".x"
1653 #else
1654         if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
1655      &    call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
1656      &      liczba(:ilen(liczba))//".cx")
1657         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
1658      &  //".cx"
1659 #endif
1660       endif
1661 #else
1662       if(mdpdb) then
1663          if (ilen(tmpdir).gt.0)
1664      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
1665          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
1666       else
1667          if (ilen(tmpdir).gt.0)
1668      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
1669          cartname=prefix(:ilen(prefix))//"_MD.cx"
1670       endif
1671 #endif
1672       if (usampl) then
1673         write (qstr,'(256(1h ))')
1674         ipos=1
1675         do i=1,nfrag
1676           iq = qinfrag(i,iset)*10
1677           iw = wfrag(i,iset)/100
1678           if (iw.gt.0) then
1679             if(me.eq.king.or..not.out1file)
1680      &       write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
1681             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
1682             ipos=ipos+7
1683           endif
1684         enddo
1685         do i=1,npair
1686           iq = qinpair(i,iset)*10
1687           iw = wpair(i,iset)/100
1688           if (iw.gt.0) then
1689             if(me.eq.king.or..not.out1file)
1690      &       write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
1691             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
1692             ipos=ipos+7
1693           endif
1694         enddo
1695 c        pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
1696 #ifdef NOXDR
1697 c        cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
1698 #else
1699 c        cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
1700 #endif
1701 c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
1702       endif
1703       icg=1
1704       if (rest) then
1705        if (restart1file) then
1706          if (me.eq.king)
1707      &     inquire(file=mremd_rst_name,exist=file_exist)
1708            write (*,*) me," Before broadcast: file_exist",file_exist
1709          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
1710      &          IERR)
1711          write (*,*) me," After broadcast: file_exist",file_exist
1712 c        inquire(file=mremd_rst_name,exist=file_exist)
1713         if(me.eq.king.or..not.out1file)
1714      &   write(iout,*) "Initial state read by master and distributed"
1715        else
1716          if (ilen(tmpdir).gt.0)
1717      &     call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
1718      &      //liczba(:ilen(liczba))//'.rst')
1719         inquire(file=rest2name,exist=file_exist)
1720        endif
1721        if(file_exist) then
1722          if(.not.restart1file) then
1723            if(me.eq.king.or..not.out1file)
1724      &      write(iout,*) "Initial state will be read from file ",
1725      &      rest2name(:ilen(rest2name))
1726            call readrst
1727          endif
1728          call rescale_weights(t_bath)
1729        else
1730         if(me.eq.king.or..not.out1file)then
1731          if (restart1file) then
1732           write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
1733      &       " does not exist"
1734          else
1735           write(iout,*) "File ",rest2name(:ilen(rest2name)),
1736      &       " does not exist"
1737          endif
1738          write(iout,*) "Initial velocities randomly generated"
1739         endif
1740         call random_vel
1741         totT=0.0d0
1742        endif
1743       else
1744 c Generate initial velocities
1745         if(me.eq.king.or..not.out1file)
1746      &   write(iout,*) "Initial velocities randomly generated"
1747         call random_vel
1748         totT=0.0d0
1749       endif
1750 c      rest2name = prefix(:ilen(prefix))//'.rst'
1751       if(me.eq.king.or..not.out1file)then
1752        write(iout,*) "Initial backbone velocities"
1753        do i=nnt,nct-1
1754         write(iout,*) (d_t(j,i),j=1,3)
1755        enddo
1756        write(iout,*) "Initial side-chain velocities"
1757        do i=nnt,nct
1758         write(iout,*) (d_t(j,i+nres),j=1,3)
1759        enddo
1760 c  Zeroing the total angular momentum of the system
1761        write(iout,*) "Calling the zero-angular
1762      & momentum subroutine"
1763       endif
1764       call inertia_tensor
1765 c  Getting the potential energy and forces and velocities and accelerations
1766       call vcm_vel(vcm)
1767 c      write (iout,*) "velocity of the center of the mass:"
1768 c      write (iout,*) (vcm(j),j=1,3)
1769       do j=1,3
1770         d_t(j,0)=d_t(j,0)-vcm(j)
1771       enddo
1772 c Removing the velocity of the center of mass
1773       call vcm_vel(vcm)
1774       if(me.eq.king.or..not.out1file)then
1775        write (iout,*) "vcm right after adjustment:"
1776        write (iout,*) (vcm(j),j=1,3)
1777       endif
1778       if (.not.rest) then
1779          call chainbuild
1780          if(iranconf.ne.0) then
1781           if (overlapsc) then
1782            print *, 'Calling OVERLAP_SC'
1783            call overlap_sc(fail)
1784           endif
1785
1786           if (searchsc) then
1787            call sc_move(2,nres-1,10,1d10,nft_sc,etot)
1788            print *,'SC_move',nft_sc,etot
1789            if(me.eq.king.or..not.out1file)
1790      &      write(iout,*) 'SC_move',nft_sc,etot
1791           endif
1792
1793           if(dccart)then
1794            print *, 'Calling MINIM_DC'
1795            call minim_dc(etot,iretcode,nfun)
1796           else
1797            call geom_to_var(nvar,varia)
1798            print *,'Calling MINIMIZE.'
1799            call minimize(etot,varia,iretcode,nfun)
1800            call var_to_geom(nvar,varia)
1801           endif
1802           if(me.eq.king.or..not.out1file)
1803      &       write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
1804          endif
1805       endif
1806       call chainbuild_cart
1807       call kinetic(EK)
1808       if (tbf) then
1809         call verlet_bath(EK)
1810       endif
1811       kinetic_T=2.0d0/(dimen*Rb)*EK
1812       if(me.eq.king.or..not.out1file)then
1813        call cartprint
1814        call intout
1815       endif
1816       call zerograd
1817       call etotal(potEcomp)
1818       potE=potEcomp(0)
1819       call cartgrad
1820       call lagrangian
1821       call max_accel
1822       if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
1823       if(me.eq.king.or..not.out1file)then
1824        write(iout,*) "Potential energy and its components"
1825        write(iout,*) (potEcomp(i),i=0,n_ene)
1826       endif
1827       potE=potEcomp(0)-potEcomp(20)
1828       totE=EK+potE
1829       itime=0
1830       if (ntwe.ne.0) call statout(itime)
1831       if(me.eq.king.or..not.out1file)
1832      &  write (iout,*) "Initial:",
1833      &   " Kinetic energy",EK," potential energy",potE,
1834      &   " total energy",totE," maximum acceleration ",
1835      &   amax
1836       if (large) then
1837         write (iout,*) "Initial coordinates"
1838         do i=1,nres
1839           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
1840      &    (c(j,i+nres),j=1,3)
1841         enddo
1842         write (iout,*) "Initial dC"
1843         do i=0,nres
1844           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
1845      &    (dc(j,i+nres),j=1,3)
1846         enddo
1847         write (iout,*) "Initial velocities"
1848         do i=0,nres
1849           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
1850      &    (d_t(j,i+nres),j=1,3)
1851         enddo
1852         write (iout,*) "Initial accelerations"
1853         do i=0,nres
1854           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
1855      &    (d_a(j,i+nres),j=1,3)
1856         enddo
1857       endif
1858       do i=0,2*nres
1859         do j=1,3
1860           dc_old(j,i)=dc(j,i)
1861           d_t_old(j,i)=d_t(j,i)
1862           d_a_old(j,i)=d_a(j,i)
1863         enddo
1864 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
1865       enddo
1866       if (RESPA) then
1867 #ifdef MPI
1868       tt0 =MPI_Wtime()
1869 #else
1870       tt0 = tcpu()
1871 #endif
1872         call zerograd
1873         call etotal_long(energia_long)
1874       if (large) write (iout,*) "energia_long",energia_long(0)
1875         call cartgrad
1876         call lagrangian
1877 #ifdef MPI
1878         t_enegrad=t_enegrad+MPI_Wtime()-tt0
1879 #else
1880         t_enegrad=t_enegrad+tcpu()-tt0
1881 #endif
1882 c        call etotal_short(energia_short)
1883 c        write (iout,*) "energia_long",energia_long(0),
1884 c     &   " energia_short",energia_short(0),
1885 c     &   " total",energia_long(0)+energia_short(0)
1886       endif
1887       return
1888       end
1889 c-----------------------------------------------------------
1890       subroutine random_vel
1891       implicit real*8 (a-h,o-z)
1892       include 'DIMENSIONS'
1893       include 'COMMON.CONTROL'
1894       include 'COMMON.VAR'
1895       include 'COMMON.MD'
1896 #ifndef LANG0
1897       include 'COMMON.LANGEVIN'
1898 #else
1899       include 'COMMON.LANGEVIN.lang0'
1900 #endif
1901       include 'COMMON.CHAIN'
1902       include 'COMMON.DERIV'
1903       include 'COMMON.GEO'
1904       include 'COMMON.LOCAL'
1905       include 'COMMON.INTERACT'
1906       include 'COMMON.IOUNITS'
1907       include 'COMMON.NAMES'
1908       include 'COMMON.TIME1'
1909       double precision xv,sigv,lowb,highb
1910 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
1911 c First generate velocities in the eigenspace of the G matrix
1912 c      write (iout,*) "Calling random_vel"
1913       xv=0.0d0
1914       do i=1,dimen
1915         sigv=dsqrt((Rb*t_bath)/geigen(i))
1916         lowb=-5*sigv
1917         highb=5*sigv
1918         d_t_work_new(i)=anorm_distr(xv,sigv,lowb,highb)
1919       enddo
1920 c      Ek1=0.0d0
1921 c      do i=1,dimen
1922 c        Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(i)**2
1923 c      enddo
1924 c Transform velocities to UNRES coordinate space
1925       do i=1,dimen
1926         d_t_work(i)=0.0d0
1927         do j=1,dimen
1928           d_t_work(i)=d_t_work(i)+Gvec(i,j)*d_t_work_new(j)
1929         enddo
1930       enddo
1931 c Transfer to the d_t vector
1932       do j=1,3
1933         d_t(j,0)=d_t_work(j)
1934       enddo
1935       ind=3
1936       do i=nnt,nct-1
1937         do j=1,3
1938           ind=ind+1
1939           if (itype(i).ne.ntyp1 .and. itype(i+1).ne.ntyp1) then
1940             d_t(j,i)=d_t_work(ind)
1941           else
1942             d_t(j,i)=0.0d0
1943           endif
1944         enddo
1945       enddo
1946       do i=nnt,nct
1947         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1948           do j=1,3
1949             ind=ind+1
1950             d_t(j,i+nres)=d_t_work(ind)
1951           enddo
1952         endif
1953       enddo
1954 c      call kinetic(EK)
1955 c      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
1956 c     &  2.0d0/(dimen*Rb)*EK,2.0d0/(dimen*Rb)*EK1
1957       return
1958       end
1959 c-----------------------------------------------------------
1960       subroutine sd_verlet_p_setup
1961 c Sets up the parameters of stochastic Verlet algorithm
1962       implicit real*8 (a-h,o-z)
1963       include 'DIMENSIONS'
1964 #ifdef MPI
1965       include 'mpif.h'
1966 #endif
1967       include 'COMMON.CONTROL'
1968       include 'COMMON.VAR'
1969       include 'COMMON.MD'
1970 #ifndef LANG0
1971       include 'COMMON.LANGEVIN'
1972 #else
1973       include 'COMMON.LANGEVIN.lang0'
1974 #endif
1975       include 'COMMON.CHAIN'
1976       include 'COMMON.DERIV'
1977       include 'COMMON.GEO'
1978       include 'COMMON.LOCAL'
1979       include 'COMMON.INTERACT'
1980       include 'COMMON.IOUNITS'
1981       include 'COMMON.NAMES'
1982       include 'COMMON.TIME1'
1983       double precision emgdt(MAXRES6),
1984      & pterm,vterm,rho,rhoc,vsig,
1985      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
1986      & afric_vec(MAXRES6),prand_vec(MAXRES6),
1987      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
1988       logical lprn /.false./
1989       double precision zero /1.0d-8/, gdt_radius /0.05d0/
1990       double precision ktm
1991 #ifdef MPI
1992       tt0 = MPI_Wtime()
1993 #else
1994       tt0 = tcpu()
1995 #endif
1996 c
1997 c AL 8/17/04 Code adapted from tinker
1998 c
1999 c Get the frictional and random terms for stochastic dynamics in the
2000 c eigenspace of mass-scaled UNRES friction matrix
2001 c
2002       do i = 1, dimen
2003             gdt = fricgam(i) * d_time
2004 c
2005 c Stochastic dynamics reduces to simple MD for zero friction
2006 c
2007             if (gdt .le. zero) then
2008                pfric_vec(i) = 1.0d0
2009                vfric_vec(i) = d_time
2010                afric_vec(i) = 0.5d0 * d_time * d_time
2011                prand_vec(i) = 0.0d0
2012                vrand_vec1(i) = 0.0d0
2013                vrand_vec2(i) = 0.0d0
2014 c
2015 c Analytical expressions when friction coefficient is large
2016 c
2017             else
2018                if (gdt .ge. gdt_radius) then
2019                   egdt = dexp(-gdt)
2020                   pfric_vec(i) = egdt
2021                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2022                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2023                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2024                   vterm = 1.0d0 - egdt**2
2025                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2026 c
2027 c Use series expansions when friction coefficient is small
2028 c
2029                else
2030                   gdt2 = gdt * gdt
2031                   gdt3 = gdt * gdt2
2032                   gdt4 = gdt2 * gdt2
2033                   gdt5 = gdt2 * gdt3
2034                   gdt6 = gdt3 * gdt3
2035                   gdt7 = gdt3 * gdt4
2036                   gdt8 = gdt4 * gdt4
2037                   gdt9 = gdt4 * gdt5
2038                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
2039      &                          - gdt5/120.0d0 + gdt6/720.0d0
2040      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
2041      &                          - gdt9/362880.0d0) / fricgam(i)**2
2042                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2043                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2044                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
2045      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
2046      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
2047      &                       + 127.0d0*gdt9/90720.0d0
2048                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
2049      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
2050      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
2051      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2052                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
2053      &                       - 17.0d0*gdt2/1280.0d0
2054      &                       + 17.0d0*gdt3/6144.0d0
2055      &                       + 40967.0d0*gdt4/34406400.0d0
2056      &                       - 57203.0d0*gdt5/275251200.0d0
2057      &                       - 1429487.0d0*gdt6/13212057600.0d0)
2058                end if
2059 c
2060 c Compute the scaling factors of random terms for the nonzero friction case
2061 c
2062                ktm = 0.5d0*d_time/fricgam(i)
2063                psig = dsqrt(ktm*pterm) / fricgam(i)
2064                vsig = dsqrt(ktm*vterm)
2065                rhoc = dsqrt(1.0d0 - rho*rho)
2066                prand_vec(i) = psig
2067                vrand_vec1(i) = vsig * rho
2068                vrand_vec2(i) = vsig * rhoc
2069             end if
2070       end do
2071       if (lprn) then
2072       write (iout,*)
2073      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2074      &  " vrand_vec2"
2075       do i=1,dimen
2076         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2077      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2078       enddo
2079       endif
2080 c
2081 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2082 c
2083 #ifndef   LANG0
2084       call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2085       call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2086       call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2087       call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2088       call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec1,vrand_mat1)
2089       call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2090 #endif
2091 c      call eigtransf1(dimen,maxres6,mt3mt2,pfric_vec,pfric_mat)
2092 c      call eigtransf1(dimen,maxres6,mt3mt2,vfric_vec,vfric_mat)
2093 c      call eigtransf1(dimen,maxres6,mt3mt2,afric_vec,afric_mat)
2094 c      call eigtransf1(dimen,maxres6,mt3mt1,prand_vec,prand_mat)
2095 c      call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec1,vrand_mat1)
2096 c      call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec2,vrand_mat2)
2097 #ifdef MPI
2098       t_sdsetup=t_sdsetup+MPI_Wtime()
2099 #else
2100       t_sdsetup=t_sdsetup+tcpu()-tt0
2101 #endif
2102       return
2103       end
2104 c-------------------------------------------------------------
2105       subroutine eigtransf1(n,ndim,ab,d,c)
2106       implicit none
2107       integer n,ndim
2108       double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
2109       integer i,j,k
2110       do i=1,n
2111         do j=1,n
2112           c(i,j)=0.0d0
2113           do k=1,n
2114             c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2115           enddo
2116         enddo
2117       enddo
2118       return
2119       end
2120 c-------------------------------------------------------------
2121       subroutine eigtransf(n,ndim,a,b,d,c)
2122       implicit none
2123       integer n,ndim
2124       double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2125       integer i,j,k
2126       do i=1,n
2127         do j=1,n
2128           c(i,j)=0.0d0
2129           do k=1,n
2130             c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2131           enddo
2132         enddo
2133       enddo
2134       return
2135       end
2136 c-------------------------------------------------------------
2137       subroutine sd_verlet1
2138 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2139       implicit real*8 (a-h,o-z)
2140       include 'DIMENSIONS'
2141       include 'COMMON.CONTROL'
2142       include 'COMMON.VAR'
2143       include 'COMMON.MD'
2144 #ifndef LANG0
2145       include 'COMMON.LANGEVIN'
2146 #else
2147       include 'COMMON.LANGEVIN.lang0'
2148 #endif
2149       include 'COMMON.CHAIN'
2150       include 'COMMON.DERIV'
2151       include 'COMMON.GEO'
2152       include 'COMMON.LOCAL'
2153       include 'COMMON.INTERACT'
2154       include 'COMMON.IOUNITS'
2155       include 'COMMON.NAMES'
2156       double precision stochforcvec(MAXRES6)
2157       common /stochcalc/ stochforcvec
2158       logical lprn /.false./
2159
2160 c      write (iout,*) "dc_old"
2161 c      do i=0,nres
2162 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2163 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2164 c      enddo
2165       do j=1,3
2166         dc_work(j)=dc_old(j,0)
2167         d_t_work(j)=d_t_old(j,0)
2168         d_a_work(j)=d_a_old(j,0)
2169       enddo
2170       ind=3
2171       do i=nnt,nct-1
2172         do j=1,3
2173           dc_work(ind+j)=dc_old(j,i)
2174           d_t_work(ind+j)=d_t_old(j,i)
2175           d_a_work(ind+j)=d_a_old(j,i)
2176         enddo
2177         ind=ind+3
2178       enddo
2179       do i=nnt,nct
2180         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2181           do j=1,3
2182             dc_work(ind+j)=dc_old(j,i+nres)
2183             d_t_work(ind+j)=d_t_old(j,i+nres)
2184             d_a_work(ind+j)=d_a_old(j,i+nres)
2185           enddo
2186           ind=ind+3
2187         endif
2188       enddo
2189 #ifndef LANG0
2190       if (lprn) then
2191       write (iout,*)
2192      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2193      &  " vrand_mat2"
2194       do i=1,dimen
2195         do j=1,dimen
2196           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2197      &      vfric_mat(i,j),afric_mat(i,j),
2198      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2199         enddo
2200       enddo
2201       endif
2202       do i=1,dimen
2203         ddt1=0.0d0
2204         ddt2=0.0d0
2205         do j=1,dimen
2206           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2207      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2208           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2209           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2210         enddo
2211         d_t_work_new(i)=ddt1+0.5d0*ddt2
2212         d_t_work(i)=ddt1+ddt2
2213       enddo
2214 #endif
2215       do j=1,3
2216         dc(j,0)=dc_work(j)
2217         d_t(j,0)=d_t_work(j)
2218       enddo
2219       ind=3
2220       do i=nnt,nct-1
2221         do j=1,3
2222           dc(j,i)=dc_work(ind+j)
2223           d_t(j,i)=d_t_work(ind+j)
2224         enddo
2225         ind=ind+3
2226       enddo
2227       do i=nnt,nct
2228         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2229           inres=i+nres
2230           do j=1,3
2231             dc(j,inres)=dc_work(ind+j)
2232             d_t(j,inres)=d_t_work(ind+j)
2233           enddo
2234           ind=ind+3
2235         endif
2236       enddo
2237       return
2238       end
2239 c--------------------------------------------------------------------------
2240       subroutine sd_verlet2
2241 c  Calculating the adjusted velocities for accelerations
2242       implicit real*8 (a-h,o-z)
2243       include 'DIMENSIONS'
2244       include 'COMMON.CONTROL'
2245       include 'COMMON.VAR'
2246       include 'COMMON.MD'
2247 #ifndef LANG0
2248       include 'COMMON.LANGEVIN'
2249 #else
2250       include 'COMMON.LANGEVIN.lang0'
2251 #endif
2252       include 'COMMON.CHAIN'
2253       include 'COMMON.DERIV'
2254       include 'COMMON.GEO'
2255       include 'COMMON.LOCAL'
2256       include 'COMMON.INTERACT'
2257       include 'COMMON.IOUNITS'
2258       include 'COMMON.NAMES'
2259       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2260       common /stochcalc/ stochforcvec
2261 c
2262 c Compute the stochastic forces which contribute to velocity change
2263 c
2264       call stochastic_force(stochforcvecV)
2265
2266 #ifndef LANG0
2267       do i=1,dimen
2268         ddt1=0.0d0
2269         ddt2=0.0d0
2270         do j=1,dimen
2271           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2272           ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
2273      &     vrand_mat2(i,j)*stochforcvecV(j)
2274         enddo
2275         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2276       enddo
2277 #endif
2278       do j=1,3
2279         d_t(j,0)=d_t_work(j)
2280       enddo
2281       ind=3
2282       do i=nnt,nct-1
2283         do j=1,3
2284           d_t(j,i)=d_t_work(ind+j)
2285         enddo
2286         ind=ind+3
2287       enddo
2288       do i=nnt,nct
2289         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2290           inres=i+nres
2291           do j=1,3
2292             d_t(j,inres)=d_t_work(ind+j)
2293           enddo
2294           ind=ind+3
2295         endif
2296       enddo
2297       return
2298       end
2299 c-----------------------------------------------------------
2300       subroutine sd_verlet_ciccotti_setup
2301 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
2302 c version
2303       implicit real*8 (a-h,o-z)
2304       include 'DIMENSIONS'
2305 #ifdef MPI
2306       include 'mpif.h'
2307 #endif
2308       include 'COMMON.CONTROL'
2309       include 'COMMON.VAR'
2310       include 'COMMON.MD'
2311 #ifndef LANG0
2312       include 'COMMON.LANGEVIN'
2313 #else
2314       include 'COMMON.LANGEVIN.lang0'
2315 #endif
2316       include 'COMMON.CHAIN'
2317       include 'COMMON.DERIV'
2318       include 'COMMON.GEO'
2319       include 'COMMON.LOCAL'
2320       include 'COMMON.INTERACT'
2321       include 'COMMON.IOUNITS'
2322       include 'COMMON.NAMES'
2323       include 'COMMON.TIME1'
2324       double precision emgdt(MAXRES6),
2325      & pterm,vterm,rho,rhoc,vsig,
2326      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
2327      & afric_vec(MAXRES6),prand_vec(MAXRES6),
2328      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
2329       logical lprn /.false./
2330       double precision zero /1.0d-8/, gdt_radius /0.05d0/
2331       double precision ktm
2332 #ifdef MPI
2333       tt0 = MPI_Wtime()
2334 #else
2335       tt0 = tcpu()
2336 #endif
2337 c
2338 c AL 8/17/04 Code adapted from tinker
2339 c
2340 c Get the frictional and random terms for stochastic dynamics in the
2341 c eigenspace of mass-scaled UNRES friction matrix
2342 c
2343       do i = 1, dimen
2344             write (iout,*) "i",i," fricgam",fricgam(i)
2345             gdt = fricgam(i) * d_time
2346 c
2347 c Stochastic dynamics reduces to simple MD for zero friction
2348 c
2349             if (gdt .le. zero) then
2350                pfric_vec(i) = 1.0d0
2351                vfric_vec(i) = d_time
2352                afric_vec(i) = 0.5d0*d_time*d_time
2353                prand_vec(i) = afric_vec(i)
2354                vrand_vec2(i) = vfric_vec(i)
2355 c
2356 c Analytical expressions when friction coefficient is large
2357 c
2358             else
2359                egdt = dexp(-gdt)
2360                pfric_vec(i) = egdt
2361                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
2362                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
2363                prand_vec(i) = afric_vec(i)
2364                vrand_vec2(i) = vfric_vec(i)
2365 c
2366 c Compute the scaling factors of random terms for the nonzero friction case
2367 c
2368 c               ktm = 0.5d0*d_time/fricgam(i)
2369 c               psig = dsqrt(ktm*pterm) / fricgam(i)
2370 c               vsig = dsqrt(ktm*vterm)
2371 c               prand_vec(i) = psig*afric_vec(i)
2372 c               vrand_vec2(i) = vsig*vfric_vec(i)
2373             end if
2374       end do
2375       if (lprn) then
2376       write (iout,*)
2377      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
2378      &  " vrand_vec2"
2379       do i=1,dimen
2380         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
2381      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2382       enddo
2383       endif
2384 c
2385 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2386 c
2387       call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
2388       call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
2389       call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
2390       call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
2391       call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
2392 #ifdef MPI
2393       t_sdsetup=t_sdsetup+MPI_Wtime()
2394 #else
2395       t_sdsetup=t_sdsetup+tcpu()-tt0
2396 #endif
2397       return
2398       end
2399 c-------------------------------------------------------------
2400       subroutine sd_verlet1_ciccotti
2401 c Applying stochastic velocity Verlet algorithm - step 1 to velocities
2402       implicit real*8 (a-h,o-z)
2403       include 'DIMENSIONS'
2404 #ifdef MPI
2405       include 'mpif.h'
2406 #endif
2407       include 'COMMON.CONTROL'
2408       include 'COMMON.VAR'
2409       include 'COMMON.MD'
2410 #ifndef LANG0
2411       include 'COMMON.LANGEVIN'
2412 #else
2413       include 'COMMON.LANGEVIN.lang0'
2414 #endif
2415       include 'COMMON.CHAIN'
2416       include 'COMMON.DERIV'
2417       include 'COMMON.GEO'
2418       include 'COMMON.LOCAL'
2419       include 'COMMON.INTERACT'
2420       include 'COMMON.IOUNITS'
2421       include 'COMMON.NAMES'
2422       double precision stochforcvec(MAXRES6)
2423       common /stochcalc/ stochforcvec
2424       logical lprn /.false./
2425
2426 c      write (iout,*) "dc_old"
2427 c      do i=0,nres
2428 c        write (iout,'(i5,3f10.5,5x,3f10.5)')
2429 c     &   i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2430 c      enddo
2431       do j=1,3
2432         dc_work(j)=dc_old(j,0)
2433         d_t_work(j)=d_t_old(j,0)
2434         d_a_work(j)=d_a_old(j,0)
2435       enddo
2436       ind=3
2437       do i=nnt,nct-1
2438         do j=1,3
2439           dc_work(ind+j)=dc_old(j,i)
2440           d_t_work(ind+j)=d_t_old(j,i)
2441           d_a_work(ind+j)=d_a_old(j,i)
2442         enddo
2443         ind=ind+3
2444       enddo
2445       do i=nnt,nct
2446         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2447           do j=1,3
2448             dc_work(ind+j)=dc_old(j,i+nres)
2449             d_t_work(ind+j)=d_t_old(j,i+nres)
2450             d_a_work(ind+j)=d_a_old(j,i+nres)
2451           enddo
2452           ind=ind+3
2453         endif
2454       enddo
2455
2456 #ifndef LANG0
2457       if (lprn) then
2458       write (iout,*)
2459      &  "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
2460      &  " vrand_mat2"
2461       do i=1,dimen
2462         do j=1,dimen
2463           write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
2464      &      vfric_mat(i,j),afric_mat(i,j),
2465      &      prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2466         enddo
2467       enddo
2468       endif
2469       do i=1,dimen
2470         ddt1=0.0d0
2471         ddt2=0.0d0
2472         do j=1,dimen
2473           dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
2474      &      +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2475           ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2476           ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2477         enddo
2478         d_t_work_new(i)=ddt1+0.5d0*ddt2
2479         d_t_work(i)=ddt1+ddt2
2480       enddo
2481 #endif
2482       do j=1,3
2483         dc(j,0)=dc_work(j)
2484         d_t(j,0)=d_t_work(j)
2485       enddo
2486       ind=3
2487       do i=nnt,nct-1
2488         do j=1,3
2489           dc(j,i)=dc_work(ind+j)
2490           d_t(j,i)=d_t_work(ind+j)
2491         enddo
2492         ind=ind+3
2493       enddo
2494       do i=nnt,nct
2495         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2496           inres=i+nres
2497           do j=1,3
2498             dc(j,inres)=dc_work(ind+j)
2499             d_t(j,inres)=d_t_work(ind+j)
2500           enddo
2501           ind=ind+3
2502         endif
2503       enddo
2504       return
2505       end
2506 c--------------------------------------------------------------------------
2507       subroutine sd_verlet2_ciccotti
2508 c  Calculating the adjusted velocities for accelerations
2509       implicit real*8 (a-h,o-z)
2510       include 'DIMENSIONS'
2511       include 'COMMON.CONTROL'
2512       include 'COMMON.VAR'
2513       include 'COMMON.MD'
2514 #ifndef LANG0
2515       include 'COMMON.LANGEVIN'
2516 #else
2517       include 'COMMON.LANGEVIN.lang0'
2518 #endif
2519       include 'COMMON.CHAIN'
2520       include 'COMMON.DERIV'
2521       include 'COMMON.GEO'
2522       include 'COMMON.LOCAL'
2523       include 'COMMON.INTERACT'
2524       include 'COMMON.IOUNITS'
2525       include 'COMMON.NAMES'
2526       double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
2527       common /stochcalc/ stochforcvec
2528 c
2529 c Compute the stochastic forces which contribute to velocity change
2530 c
2531       call stochastic_force(stochforcvecV)
2532 #ifndef LANG0
2533       do i=1,dimen
2534         ddt1=0.0d0
2535         ddt2=0.0d0
2536         do j=1,dimen
2537
2538           ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2539 c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
2540           ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
2541         enddo
2542         d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2543       enddo
2544 #endif
2545       do j=1,3
2546         d_t(j,0)=d_t_work(j)
2547       enddo
2548       ind=3
2549       do i=nnt,nct-1
2550         do j=1,3
2551           d_t(j,i)=d_t_work(ind+j)
2552         enddo
2553         ind=ind+3
2554       enddo
2555       do i=nnt,nct
2556         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2557           inres=i+nres
2558           do j=1,3
2559             d_t(j,inres)=d_t_work(ind+j)
2560           enddo
2561           ind=ind+3
2562         endif
2563       enddo
2564       return
2565       end