source/unres/src_MD/src/md-diff/np/a

   1 c------------------------------------------------------
   2       double precision function HNose(ek,s,e,pi,Q,t_bath,dimen)
   3       implicit none
   4       double precision ek,s,e,pi,Q,t_bath,Rb
   5       integer dimen
   6       Rb=0.001986d0
   7       HNose=ek+e+pi**2/(2*Q)+dimen*Rb*t_bath*log(s)
   8 c      print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimen,"--",
   9 c     &      pi**2/(2*Q),dimen*Rb*t_bath*log(s)
  10       return
  11       end
  12 c-----------------------------------------------------------------
  13       double precision function HNose_nh(eki,e)
  14       implicit real*8 (a-h,o-z)
  15       include 'DIMENSIONS'
  16       include 'COMMON.MD'
  17       HNose_nh=eki+e+dimen*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2
  18       do i=2,nnos
  19         HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i)
  20       enddo
  21 c      write(4,'(5e15.5)')
  22 c     &       vlogs(1),xlogs(1),HNose,eki,e
  23       return
  24       end
  25 c-----------------------------------------------------------------
  26       SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
  27       implicit real*8 (a-h,o-z)
  28       include 'DIMENSIONS'
  29       include 'COMMON.MD'
  30       double precision akin,gnkt,dt,aa,gkt,scale
  31       double precision wdti(maxyosh),wdti2(maxyosh),
  32      &                 wdti4(maxyosh),wdti8(maxyosh)
  33       integer i,iresn,iyosh,inos,nnos1
  34
  35       dt=d_time
  36       nnos1=nnos+1
  37       GKT = Rb*t_bath
  38       GNKT = dimen*GKT
  39       akin=akin*2
  40
  41
  42 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
  43 C INTEGRATION FROM t=0 TO t=DT/2
  44 C GET THE TOTAL KINETIC ENERGY
  45       SCALE = 1.D0
  46 c      CALL GETKINP(MASS,VX,VY,VZ,AKIN)
  47 C UPDATE THE FORCES
  48       GLOGS(1) = (AKIN - GNKT)/QMASS(1)
  49 C START THE MULTIPLE TIME STEP PROCEDURE
  50       DO IRESN = 1,NRESN
  51        DO IYOSH = 1,NYOSH
  52 C UPDATE THE THERMOSTAT VELOCITIES
  53         VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
  54         DO INOS = 1,NNOS-1
  55          AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
  56          VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
  57      &          + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
  58         ENDDO
  59 C UPDATE THE PARTICLE VELOCITIES
  60         AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
  61         SCALE = SCALE*AA
  62 C UPDATE THE FORCES
  63         GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
  64 C UPDATE THE THERMOSTAT POSITIONS
  65         DO INOS = 1,NNOS
  66          XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
  67         ENDDO
  68 C UPDATE THE THERMOSTAT VELOCITIES
  69         DO INOS = 1,NNOS-1
  70          AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
  71          VLOGS(INOS) = VLOGS(INOS)*AA*AA
  72      &      + WDTI4(IYOSH)*GLOGS(INOS)*AA
  73          GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
  74      &      -GKT)/QMASS(INOS+1)
  75         ENDDO
  76         VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
  77        ENDDO
  78       ENDDO
  79 C UPDATE THE PARTICLE VELOCITIES
  80 c outside of this subroutine
  81 c      DO I = 1,N
  82 c       VX(I) = VX(I)*SCALE
  83 c       VY(I) = VY(I)*SCALE
  84 c       VZ(I) = VZ(I)*SCALE
  85 c      ENDDO
  86       RETURN
  87       END
  88 c-----------------------------------------------------------------
  89       subroutine tnp1_respa_i_step1
  90 c Applying Nose-Poincare algorithm - step 1 to coordinates
  91 c JPSJ 70 75 (2001) S. Nose
  92 c
  93 c d_t is not updated here
  94 c
  95       implicit real*8 (a-h,o-z)
  96       include 'DIMENSIONS'
  97       include 'COMMON.CONTROL'
  98       include 'COMMON.VAR'
  99       include 'COMMON.MD'
 100       include 'COMMON.CHAIN'
 101       include 'COMMON.DERIV'
 102       include 'COMMON.GEO'
 103       include 'COMMON.LOCAL'
 104       include 'COMMON.INTERACT'
 105       include 'COMMON.IOUNITS'
 106       include 'COMMON.NAMES'
 107       double precision adt,adt2,tmp
 108
 109       tmp=1+pi_np/(2*Q_np)*0.5*d_time
 110       s12_np=s_np*tmp**2
 111       pistar=pi_np/tmp
 112       s12_dt=d_time/s12_np
 113       d_time_s12=d_time*0.5*s12_np
 114
 115       do j=1,3
 116         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
 117         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
 118       enddo
 119       do i=nnt,nct-1
 120         do j=1,3
 121           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
 122           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
 123         enddo
 124       enddo
 125       do i=nnt,nct
 126         if (itype(i).ne.10) then
 127           inres=i+nres
 128           do j=1,3
 129            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
 130            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
 131           enddo
 132         endif
 133       enddo
 134       return
 135       end
 136 c---------------------------------------------------------------------
 137       subroutine tnp1_respa_i_step2
 138 c  Step 2 of the velocity Verlet algorithm: update velocities
 139       implicit real*8 (a-h,o-z)
 140       include 'DIMENSIONS'
 141       include 'COMMON.CONTROL'
 142       include 'COMMON.VAR'
 143       include 'COMMON.MD'
 144       include 'COMMON.CHAIN'
 145       include 'COMMON.DERIV'
 146       include 'COMMON.GEO'
 147       include 'COMMON.LOCAL'
 148       include 'COMMON.INTERACT'
 149       include 'COMMON.IOUNITS'
 150       include 'COMMON.NAMES'
 151
 152       double precision d_time_s12
 153
 154       do i=0,2*nres
 155        do j=1,3
 156         d_t(j,i)=d_t_new(j,i)
 157        enddo
 158       enddo
 159
 160       call kinetic(EK)
 161       EK=EK/s12_np**2
 162
 163       d_time_s12=0.5d0*s12_np*d_time
 164
 165       do j=1,3
 166         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
 167       enddo
 168       do i=nnt,nct-1
 169         do j=1,3
 170           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
 171         enddo
 172       enddo
 173       do i=nnt,nct
 174         if (itype(i).ne.10) then
 175           inres=i+nres
 176           do j=1,3
 177             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
 178           enddo
 179         endif
 180       enddo
 181
 182       pistar=pistar+(EK-0.5*(E_old+potE)
 183      &       -dimen*Rb*t_bath*log(s12_np)+Csplit-dimen*Rb*t_bath)*d_time
 184       tmp=1+pistar/(2*Q_np)*0.5*d_time
 185       s_np=s12_np*tmp**2
 186       pi_np=pistar/tmp
 187
 188       return
 189       end
 190 c-------------------------------------------------------
 191
 192       subroutine tnp1_step1
 193 c Applying Nose-Poincare algorithm - step 1 to coordinates
 194 c JPSJ 70 75 (2001) S. Nose
 195 c
 196 c d_t is not updated here
 197 c
 198       implicit real*8 (a-h,o-z)
 199       include 'DIMENSIONS'
 200       include 'COMMON.CONTROL'
 201       include 'COMMON.VAR'
 202       include 'COMMON.MD'
 203       include 'COMMON.CHAIN'
 204       include 'COMMON.DERIV'
 205       include 'COMMON.GEO'
 206       include 'COMMON.LOCAL'
 207       include 'COMMON.INTERACT'
 208       include 'COMMON.IOUNITS'
 209       include 'COMMON.NAMES'
 210       double precision adt,adt2,tmp
 211
 212       tmp=1+pi_np/(2*Q_np)*0.5*d_time
 213       s12_np=s_np*tmp**2
 214       pistar=pi_np/tmp
 215       s12_dt=d_time/s12_np
 216       d_time_s12=d_time*0.5*s12_np
 217
 218       do j=1,3
 219         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
 220         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
 221       enddo
 222       do i=nnt,nct-1
 223         do j=1,3
 224           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
 225           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
 226         enddo
 227       enddo
 228       do i=nnt,nct
 229         if (itype(i).ne.10) then
 230           inres=i+nres
 231           do j=1,3
 232            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
 233            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
 234           enddo
 235         endif
 236       enddo
 237       return
 238       end
 239 c---------------------------------------------------------------------
 240       subroutine tnp1_step2
 241 c  Step 2 of the velocity Verlet algorithm: update velocities
 242       implicit real*8 (a-h,o-z)
 243       include 'DIMENSIONS'
 244       include 'COMMON.CONTROL'
 245       include 'COMMON.VAR'
 246       include 'COMMON.MD'
 247       include 'COMMON.CHAIN'
 248       include 'COMMON.DERIV'
 249       include 'COMMON.GEO'
 250       include 'COMMON.LOCAL'
 251       include 'COMMON.INTERACT'
 252       include 'COMMON.IOUNITS'
 253       include 'COMMON.NAMES'
 254
 255       double precision d_time_s12
 256
 257       do i=0,2*nres
 258        do j=1,3
 259         d_t(j,i)=d_t_new(j,i)
 260        enddo
 261       enddo
 262
 263       call kinetic(EK)
 264       EK=EK/s12_np**2
 265
 266       d_time_s12=0.5d0*s12_np*d_time
 267
 268       do j=1,3
 269         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
 270       enddo
 271       do i=nnt,nct-1
 272         do j=1,3
 273           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
 274         enddo
 275       enddo
 276       do i=nnt,nct
 277         if (itype(i).ne.10) then
 278           inres=i+nres
 279           do j=1,3
 280             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
 281           enddo
 282         endif
 283       enddo
 284
 285 cd      write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
 286       pistar=pistar+(EK-0.5*(E_old+potE)
 287      &       -dimen*Rb*t_bath*log(s12_np)+H0-dimen*Rb*t_bath)*d_time
 288       tmp=1+pistar/(2*Q_np)*0.5*d_time
 289       s_np=s12_np*tmp**2
 290       pi_np=pistar/tmp
 291
 292       return
 293       end
 294
 295 c-----------------------------------------------------------------
 296       subroutine tnp_respa_i_step1
 297 c Applying Nose-Poincare algorithm - step 1 to coordinates
 298 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
 299 c
 300 c d_t is not updated here, it is destroyed
 301 c
 302       implicit real*8 (a-h,o-z)
 303       include 'DIMENSIONS'
 304       include 'COMMON.CONTROL'
 305       include 'COMMON.VAR'
 306       include 'COMMON.MD'
 307       include 'COMMON.CHAIN'
 308       include 'COMMON.DERIV'
 309       include 'COMMON.GEO'
 310       include 'COMMON.LOCAL'
 311       include 'COMMON.INTERACT'
 312       include 'COMMON.IOUNITS'
 313       include 'COMMON.NAMES'
 314       double precision C_np,d_time_s,tmp,d_time_ss
 315
 316       d_time_s=d_time*0.5*s_np
 317 ct2      d_time_s=d_time*0.5*s12_np
 318
 319       do j=1,3
 320         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
 321       enddo
 322       do i=nnt,nct-1
 323         do j=1,3
 324           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
 325         enddo
 326       enddo
 327       do i=nnt,nct
 328         if (itype(i).ne.10) then
 329           inres=i+nres
 330           do j=1,3
 331            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
 332           enddo
 333         endif
 334       enddo
 335
 336       do i=0,2*nres
 337        do j=1,3
 338         d_t(j,i)=d_t_new(j,i)
 339        enddo
 340       enddo
 341
 342       call kinetic(EK)
 343       EK=EK/s_np**2
 344
 345       C_np=0.5*d_time*(dimen*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit)
 346      &                     -pi_np
 347
 348       pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
 349       tmp=0.5*d_time*pistar/Q_np
 350       s12_np=s_np*(1.0+tmp)/(1.0-tmp)
 351
 352       d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
 353 ct2      d_time_ss=d_time/s12_np
 354 c      d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np)
 355
 356       do j=1,3
 357         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
 358       enddo
 359       do i=nnt,nct-1
 360         do j=1,3
 361           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
 362         enddo
 363       enddo
 364       do i=nnt,nct
 365         if (itype(i).ne.10) then
 366           inres=i+nres
 367           do j=1,3
 368            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
 369           enddo
 370         endif
 371       enddo
 372
 373       return
 374       end
 375 c---------------------------------------------------------------------
 376
 377       subroutine tnp_respa_i_step2
 378 c  Step 2 of the velocity Verlet algorithm: update velocities
 379       implicit real*8 (a-h,o-z)
 380       include 'DIMENSIONS'
 381       include 'COMMON.CONTROL'
 382       include 'COMMON.VAR'
 383       include 'COMMON.MD'
 384       include 'COMMON.CHAIN'
 385       include 'COMMON.DERIV'
 386       include 'COMMON.GEO'
 387       include 'COMMON.LOCAL'
 388       include 'COMMON.INTERACT'
 389       include 'COMMON.IOUNITS'
 390       include 'COMMON.NAMES'
 391
 392       double precision d_time_s
 393
 394       EK=EK*(s_np/s12_np)**2
 395       HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen)
 396       pi_np=pistar+0.5*d_time*(2*EK-dimen*Rb*t_bath
 397      &                              -HNose1+Csplit)
 398
 399 cr      print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long
 400       d_time_s=d_time*0.5*s12_np
 401 c      d_time_s=d_time*0.5*s_np
 402
 403       do j=1,3
 404         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
 405       enddo
 406       do i=nnt,nct-1
 407         do j=1,3
 408           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
 409         enddo
 410       enddo
 411       do i=nnt,nct
 412         if (itype(i).ne.10) then
 413           inres=i+nres
 414           do j=1,3
 415             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
 416           enddo
 417         endif
 418       enddo
 419
 420       s_np=s12_np
 421
 422       return
 423       end
 424 c-----------------------------------------------------------------
 425       subroutine tnp_respa_step1
 426 c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
 427 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
 428 c
 429 c d_t is not updated here, it is destroyed
 430 c
 431       implicit real*8 (a-h,o-z)
 432       include 'DIMENSIONS'
 433       include 'COMMON.CONTROL'
 434       include 'COMMON.VAR'
 435       include 'COMMON.MD'
 436       include 'COMMON.CHAIN'
 437       include 'COMMON.DERIV'
 438       include 'COMMON.GEO'
 439       include 'COMMON.LOCAL'
 440       include 'COMMON.INTERACT'
 441       include 'COMMON.IOUNITS'
 442       include 'COMMON.NAMES'
 443       double precision C_np,d_time_s,tmp,d_time_ss
 444       double precision energia(0:n_ene)
 445
 446       d_time_s=d_time*0.5*s_np
 447
 448       do j=1,3
 449         d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
 450       enddo
 451       do i=nnt,nct-1
 452         do j=1,3
 453           d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
 454         enddo
 455       enddo
 456       do i=nnt,nct
 457         if (itype(i).ne.10) then
 458           inres=i+nres
 459           do j=1,3
 460            d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
 461           enddo
 462         endif
 463       enddo
 464
 465
 466 c      C_np=0.5*d_time*(dimen*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
 467 c     &                     -pi_np
 468 c
 469 c      pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
 470 c      tmp=0.5*d_time*pistar/Q_np
 471 c      s12_np=s_np*(1.0+tmp)/(1.0-tmp)
 472 c      write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
 473
 474 ct1      pi_np=pistar
 475 c      sold_np=s_np
 476 c      s_np=s12_np
 477
 478 c-------------------------------------
 479 c test of reviewer's comment
 480        pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
 481 cr       print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long
 482 c-------------------------------------
 483
 484       return
 485       end
 486 c---------------------------------------------------------------------
 487       subroutine tnp_respa_step2
 488 c  Step 2 of the velocity Verlet algorithm: update velocities for RESPA
 489       implicit real*8 (a-h,o-z)
 490       include 'DIMENSIONS'
 491       include 'COMMON.CONTROL'
 492       include 'COMMON.VAR'
 493       include 'COMMON.MD'
 494       include 'COMMON.CHAIN'
 495       include 'COMMON.DERIV'
 496       include 'COMMON.GEO'
 497       include 'COMMON.LOCAL'
 498       include 'COMMON.INTERACT'
 499       include 'COMMON.IOUNITS'
 500       include 'COMMON.NAMES'
 501
 502       double precision d_time_s
 503
 504 ct1      s12_np=s_np
 505 ct2      pistar=pi_np
 506
 507 ct      call kinetic(EK)
 508 ct      HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen)
 509 ct      pi_np=pistar+0.5*d_time*(2*EK-dimen*Rb*t_bath)
 510 ct     &                              -0.5*d_time*(HNose1-H0)
 511
 512 c-------------------------------------
 513 c test of reviewer's comment
 514       pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
 515 cr      print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long
 516 c-------------------------------------
 517       d_time_s=d_time*0.5*s_np
 518
 519       do j=1,3
 520         d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
 521       enddo
 522       do i=nnt,nct-1
 523         do j=1,3
 524           d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
 525         enddo
 526       enddo
 527       do i=nnt,nct
 528         if (itype(i).ne.10) then
 529           inres=i+nres
 530           do j=1,3
 531             d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
 532           enddo
 533         endif
 534       enddo
 535
 536 cd      s_np=s12_np
 537
 538       return
 539       end
 540 c---------------------------------------------------------------------
 541       subroutine tnp_step1
 542 c Applying Nose-Poincare algorithm - step 1 to coordinates
 543 c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
 544 c
 545 c d_t is not updated here, it is destroyed
 546 c
 547       implicit real*8 (a-h,o-z)
 548       include 'DIMENSIONS'
 549       include 'COMMON.CONTROL'
 550       include 'COMMON.VAR'
 551       include 'COMMON.MD'
 552       include 'COMMON.CHAIN'
 553       include 'COMMON.DERIV'
 554       include 'COMMON.GEO'
 555       include 'COMMON.LOCAL'
 556       include 'COMMON.INTERACT'
 557       include 'COMMON.IOUNITS'
 558       include 'COMMON.NAMES'
 559       double precision C_np,d_time_s,tmp,d_time_ss
 560
 561       d_time_s=d_time*0.5*s_np
 562
 563       do j=1,3
 564         d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
 565       enddo
 566       do i=nnt,nct-1
 567         do j=1,3
 568           d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
 569         enddo
 570       enddo
 571       do i=nnt,nct
 572         if (itype(i).ne.10) then
 573           inres=i+nres
 574           do j=1,3
 575            d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
 576           enddo
 577         endif
 578       enddo
 579
 580       do i=0,2*nres
 581        do j=1,3
 582         d_t(j,i)=d_t_new(j,i)
 583        enddo
 584       enddo
 585
 586       call kinetic(EK)
 587       EK=EK/s_np**2
 588
 589       C_np=0.5*d_time*(dimen*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
 590      &                     -pi_np
 591
 592       pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
 593       tmp=0.5*d_time*pistar/Q_np
 594       s12_np=s_np*(1.0+tmp)/(1.0-tmp)
 595 c      write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
 596
 597       d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
 598
 599       do j=1,3
 600         dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
 601       enddo
 602       do i=nnt,nct-1
 603         do j=1,3
 604           dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
 605         enddo
 606       enddo
 607       do i=nnt,nct
 608         if (itype(i).ne.10) then
 609           inres=i+nres
 610           do j=1,3
 611            dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
 612           enddo
 613         endif
 614       enddo
 615
 616       return
 617       end
 618 c-----------------------------------------------------------------
 619       subroutine tnp_step2
 620 c  Step 2 of the velocity Verlet algorithm: update velocities
 621       implicit real*8 (a-h,o-z)
 622       include 'DIMENSIONS'
 623       include 'COMMON.CONTROL'
 624       include 'COMMON.VAR'
 625       include 'COMMON.MD'
 626       include 'COMMON.CHAIN'
 627       include 'COMMON.DERIV'
 628       include 'COMMON.GEO'
 629       include 'COMMON.LOCAL'
 630       include 'COMMON.INTERACT'
 631       include 'COMMON.IOUNITS'
 632       include 'COMMON.NAMES'
 633
 634       double precision d_time_s
 635
 636       EK=EK*(s_np/s12_np)**2
 637       HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen)
 638       pi_np=pistar+0.5*d_time*(2*EK-dimen*Rb*t_bath)
 639      &                              -0.5*d_time*(HNose1-H0)
 640
 641 cd      write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
 642       d_time_s=d_time*0.5*s12_np
 643
 644       do j=1,3
 645         d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
 646       enddo
 647       do i=nnt,nct-1
 648         do j=1,3
 649           d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
 650         enddo
 651       enddo
 652       do i=nnt,nct
 653         if (itype(i).ne.10) then
 654           inres=i+nres
 655           do j=1,3
 656             d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
 657           enddo
 658         endif
 659       enddo
 660
 661       s_np=s12_np
 662
 663       return
 664       end