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