1* 2* $Id$ 3* 4 5* *************************** 6* * * 7* * Nose_Init * 8* * * 9* *************************** 10 11 subroutine Nose_Init(nemax,eke0) 12 implicit none 13 integer nemax 14 real*8 eke0 15 16 17#include "bafdecls.fh" 18#include "btdb.fh" 19#include "errquit.fh" 20#include "nose-hoover.fh" 21 22 23* **** boltzman constant **** 24 double precision kb 25 parameter (kb=3.16679d-6) 26 27 28* **** local variables **** 29 logical value,nosers 30 integer i,rtdb,n,m 31 real*8 pi,am,fmass,tmpe 32 real*8 Te,Pe_tmp,betae 33 real*8 Tr,Pr_tmp 34 35* **** external functions **** 36 real*8 control_fake_mass 37 real*8 control_Nose_Pe 38 real*8 control_Nose_Te 39 real*8 control_Nose_Pr 40 real*8 control_Nose_Tr 41 real*8 ion_amass 42 integer ion_nion,control_rtdb 43 44 external control_fake_mass 45 external control_Nose_Pe 46 external control_Nose_Te 47 external control_Nose_Pr 48 external control_Nose_Tr 49 external ion_amass 50 external ion_nion,control_rtdb 51 52 53* ************************************ 54* **** initialize the thermostats **** 55* ************************************ 56 Te = control_Nose_Te() 57 Tr = control_Nose_Tr() 58 Pe_tmp = control_Nose_Pe() 59 Pr_tmp = control_Nose_Pr() 60 61 rtdb = control_rtdb() 62 if (.not.btdb_get(rtdb,'cpmd:nose_restart',mt_log,1,nosers)) 63 > nosers = .true. 64 65 if (.not.btdb_get(rtdb,'cpmd:Mchain',mt_int,1,Mchain)) Mchain = 1 66 if (.not.btdb_get(rtdb,'cpmd:Nchain',mt_int,1,Nchain)) Nchain = 1 67 if (.not.btdb_get(rtdb,'cpmd:Ne_chain',mt_dbl,1,Ne_chain)) 68 > Ne_chain = 3.0d0*nemax 69 if ( (.not.btdb_get(rtdb,'cpmd:eke0',mt_dbl,1,eke0_init)) 70 > .or.(.not.nosers)) 71 > eke0_init = eke0 72 73 fmass = control_fake_mass() 74 75 value = BA_alloc_get(mt_dbl,Mchain,'Xem',Xem(2),Xem(1)) 76 value = value.and.BA_alloc_get(mt_dbl,Mchain,'Xe0',Xe0(2),Xe0(1)) 77 value = value.and.BA_alloc_get(mt_dbl,Mchain,'Xe1',Xe1(2),Xe1(1)) 78 value = value.and.BA_alloc_get(mt_dbl,Mchain,'Xe2',Xe2(2),Xe2(1)) 79 value = value.and.BA_alloc_get(mt_dbl,Mchain,'Ee0',Ee0(2),Ee0(1)) 80 value = value.and.BA_alloc_get(mt_dbl,Mchain,'Qe',Qe(2),Qe(1)) 81 value = value.and.BA_alloc_get(mt_dbl,Mchain,'Pe',Pe(2),Pe(1)) 82 83 value = value.and.BA_alloc_get(mt_dbl,Nchain,'Xrm',Xrm(2),Xrm(1)) 84 value = value.and.BA_alloc_get(mt_dbl,Nchain,'Xr0',Xr0(2),Xr0(1)) 85 value = value.and.BA_alloc_get(mt_dbl,Nchain,'Xr1',Xr1(2),Xr1(1)) 86 value = value.and.BA_alloc_get(mt_dbl,Nchain,'Xr2',Xr2(2),Xr2(1)) 87 value = value.and.BA_alloc_get(mt_dbl,Nchain,'Er0',Er0(2),Er0(1)) 88 value = value.and.BA_alloc_get(mt_dbl,Nchain,'Qr',Qr(2),Qr(1)) 89 value = value.and.BA_alloc_get(mt_dbl,Nchain,'Pr',Pr(2),Pr(1)) 90 if (.not.value) 91 > call errquit("Nose_Init: out of heap memory",0,MA_ERR) 92 93* **** restart using Newton Step ??? **** 94c call dcopy(Mchain,0.0d0,0,dbl_mb(Xem(1)),1) 95c call dcopy(Mchain,0.0d0,0,dbl_mb(Xe0(1)),1) 96c call dcopy(Mchain,0.0d0,0,dbl_mb(Xe1(1)),1) 97 if ((.not.btdb_get(rtdb,'cpmd:Xe1',mt_dbl,Mchain,dbl_mb(Xe1(1)))) 98 > .or.(.not.nosers)) 99 > call dcopy(Mchain,0.0d0,0,dbl_mb(Xe1(1)),1) 100 if ((.not.btdb_get(rtdb,'cpmd:Xe0',mt_dbl,Mchain,dbl_mb(Xe0(1)))) 101 > .or.(.not.nosers)) 102 > call dcopy(Mchain,0.0d0,0,dbl_mb(Xe0(1)),1) 103 if ((.not.btdb_get(rtdb,'cpmd:Xem',mt_dbl,Mchain,dbl_mb(Xem(1)))) 104 > .or.(.not.nosers)) 105 > call dcopy(Mchain,0.0d0,0,dbl_mb(Xem(1)),1) 106 107 call dcopy(Mchain,0.0d0,0,dbl_mb(Xe2(1)),1) 108 call dcopy(Mchain,Pe_tmp,0,dbl_mb(Pe(1)),1) 109 110 111c call dcopy(Nchain,0.0d0,0,dbl_mb(Xrm(1)),1) 112c call dcopy(Nchain,0.0d0,0,dbl_mb(Xr0(1)),1) 113c call dcopy(Nchain,0.0d0,0,dbl_mb(Xr1(1)),1) 114 if ((.not.btdb_get(rtdb,'cpmd:Xr1',mt_dbl,Nchain,dbl_mb(Xr1(1)))) 115 > .or.(.not.nosers)) 116 > call dcopy(Nchain,0.0d0,0,dbl_mb(Xr1(1)),1) 117 if ((.not.btdb_get(rtdb,'cpmd:Xr0',mt_dbl,Nchain,dbl_mb(Xr0(1)))) 118 > .or.(.not.nosers)) 119 > call dcopy(Nchain,0.0d0,0,dbl_mb(Xr0(1)),1) 120 if ((.not.btdb_get(rtdb,'cpmd:Xrm',mt_dbl,Nchain,dbl_mb(Xrm(1)))) 121 > .or.(.not.nosers)) 122 > call dcopy(Nchain,0.0d0,0,dbl_mb(Xrm(1)),1) 123 124 call dcopy(Nchain,0.0d0,0,dbl_mb(Xr2(1)),1) 125 call dcopy(Nchain,Pr_tmp,0,dbl_mb(Pr(1)),1) 126 127 128 129c Xe0 = 0.0d0 130c Xe1 = 0.0d0 131c Xe2 = 0.0d0 132c Xr0 = 0.0d0 133c Xr1 = 0.0d0 134c Xr2 = 0.0d0 135 136 137* **** Set Er0(1) = (1/2)*(g*k*T), where g=number of degrees of freedom **** 138 if (ion_nion().gt.2) then 139c Er0 = 0.5d0*(3.0d0*dble(ion_nion())-6.0d0)*kb*Tr 140 dbl_mb(Er0(1)) = 0.5d0*(3.0d0*dble(ion_nion())-6.0d0)*kb*Tr 141 142* **** Dimer molecule. note that the above formula may not work for **** 143* **** linear molecules with more than two atoms. **** 144 else 145c Er0 = 0.5d0*(1)*kb*Tr 146 dbl_mb(Er0(1)) = 0.5d0*(1)*kb*Tr 147 end if 148 149* **** Set Er0(2:Nchain) = 1/2*(k*T) **** 150 if (Nchain.gt.1) then 151 call dcopy(Nchain-1,0.5d0*kb*Tr,0,dbl_mb(Er0(1)+1),1) 152 end if 153 154 155* *** total mass *** 156 am = 0.0d0 157 do i=1,ion_nion() 158 am = am + ion_amass(i) 159 end do 160 161 162* **** Set Ee0(1) - read from rtdb otherwise use current KE **** 163 dbl_mb(Ee0(1)) = 4.0d0*kb*Te*fmass*dble(ion_nion())/am * eke0_init 164 165* **** Set Ee0(2:Mchain) = 1/2*(1/betae), where 1/betae = 2*Ee/Ne **** 166 if (Mchain.gt.1) then 167 betae = dbl_mb(Ee0(1))/Ne_chain 168 call dcopy(Mchain-1,betae,0,dbl_mb(Ee0(1)+1),1) 169 end if 170 171 172* **** Set Qe and Qr - read from rtdb otherwise set using periods **** 173 pi = 4.0d0*datan(1.0d0) 174 value = btdb_get(rtdb,'cpmd:Qe',mt_dbl,Mchain,dbl_mb(Qe(1))) 175 > .and.btdb_get(rtdb,'cpmd:Qr',mt_dbl,Nchain,dbl_mb(Qr(1))) 176 177 if ((.not.value).or.(.not.nosers)) then 178 do m=1,Mchain 179 dbl_mb(Qe(1)+m-1)=dbl_mb(Ee0(1)+m-1)*(dbl_mb(Pe(1)+m-1)/pi)**2 180 end do 181 do n=1,Nchain 182 dbl_mb(Qr(1)+n-1)=dbl_mb(Er0(1)+n-1)*(dbl_mb(Pr(1)+n-1)/pi)**2 183 end do 184 else 185 do m=1,Mchain 186 dbl_mb(Pe(1)+m-1) 187 > = pi*dsqrt(dbl_mb(Qe(1)+m-1)/dbl_mb(Ee0(1)+m-1)) 188 end do 189 do n=1,Nchain 190 dbl_mb(Pr(1)+N-1) 191 > = pi*dsqrt(dbl_mb(Qr(1)+n-1)/dbl_mb(Er0(1)+n-1)) 192 end do 193 end if 194 195 return 196 end 197 198 199* *************************** 200* * * 201* * Nose_end * 202* * * 203* *************************** 204 205 subroutine Nose_end() 206 implicit none 207 208#include "bafdecls.fh" 209#include "btdb.fh" 210#include "errquit.fh" 211#include "nose-hoover.fh" 212 213 214 !**** local variables **** 215 integer MASTER,taskid 216 parameter (MASTER=0) 217 logical value,value2 218 integer rtdb,m,n 219 real*8 dt 220 221 !**** external functions **** 222 integer control_rtdb 223 real*8 control_time_step 224 external control_rtdb 225 external control_time_step 226 227* **** put velecities in Xem and Xrm **** 228 dt = control_time_step() 229 do m=1,Mchain 230 dbl_mb(Xem(1)+m-1) = (dbl_mb(Xe2(1)+m-1)-dbl_mb(Xe0(1)+m-1)) 231 > /(2.0d0*dt) 232 end do 233 do n=1,Nchain 234 dbl_mb(Xrm(1)+n-1) = (dbl_mb(Xr2(1)+n-1)-dbl_mb(Xr0(1)+n-1)) 235 > /(2.0d0*dt) 236 end do 237 238* **** save restart information to rtdb **** 239 call Parallel_taskid(taskid) 240 rtdb = control_rtdb() 241 value2 = btdb_parallel(.false.) 242 if (taskid.eq.MASTER) then 243 value = btdb_put(rtdb,'cpmd:eke0',mt_dbl,1,eke0_init) 244 > .and.btdb_put(rtdb,'cpmd:Mchain',mt_int,1,Mchain) 245 > .and.btdb_put(rtdb,'cpmd:Nchain',mt_int,1,Nchain) 246 > .and.btdb_put(rtdb,'cpmd:Ne_chain',mt_dbl,1,Ne_chain) 247 > .and.btdb_put(rtdb,'cpmd:Qe', mt_dbl,Mchain,dbl_mb(Qe(1))) 248 > .and.btdb_put(rtdb,'cpmd:Xe1',mt_dbl,Mchain,dbl_mb(Xe1(1))) 249 > .and.btdb_put(rtdb,'cpmd:Xe0',mt_dbl,Mchain,dbl_mb(Xe0(1))) 250 > .and.btdb_put(rtdb,'cpmd:Xem',mt_dbl,Mchain,dbl_mb(Xem(1))) 251 > .and.btdb_put(rtdb,'cpmd:Qr', mt_dbl,Nchain,dbl_mb(Qr(1))) 252 > .and.btdb_put(rtdb,'cpmd:Xr1',mt_dbl,Nchain,dbl_mb(Xr1(1))) 253 > .and.btdb_put(rtdb,'cpmd:Xr0',mt_dbl,Nchain,dbl_mb(Xr0(1))) 254 > .and.btdb_put(rtdb,'cpmd:Xrm',mt_dbl,Nchain,dbl_mb(Xrm(1))) 255 end if 256 value2 = btdb_parallel(.true.) 257 if (.not.value) 258 > call errquit( 259 > 'Nose_End: error writing to rtdb', 0, RTDB_ERR) 260 261 262 value = value.and.BA_free_heap(Xem(2)) 263 value = value.and.BA_free_heap(Xe0(2)) 264 value = value.and.BA_free_heap(Xe1(2)) 265 value = value.and.BA_free_heap(Xe2(2)) 266 value = value.and.BA_free_heap(Ee0(2)) 267 value = value.and.BA_free_heap(Qe(2)) 268 value = value.and.BA_free_heap(Pe(2)) 269 270 value = value.and.BA_free_heap(Xrm(2)) 271 value = value.and.BA_free_heap(Xr0(2)) 272 value = value.and.BA_free_heap(Xr1(2)) 273 value = value.and.BA_free_heap(Xr2(2)) 274 value = value.and.BA_free_heap(Er0(2)) 275 value = value.and.BA_free_heap(Qr(2)) 276 value = value.and.BA_free_heap(Pr(2)) 277 if (.not.value) 278 > call errquit('Nose_End: error freeing heap',0,MA_ERR) 279 280 return 281 end 282 283 284 285* *************************** 286* * * 287* * Nose_reset_T * 288* * * 289* *************************** 290 291 subroutine Nose_reset_T(Te_new,Tr_new) 292 implicit none 293 real*8 Te_new,Tr_new 294 295#include "bafdecls.fh" 296#include "nose-hoover.fh" 297 298* **** boltzman constant **** 299 double precision kb 300 parameter (kb=3.16679d-6) 301 302* **** local variables **** 303 integer i 304 real*8 am,fmass,betae 305 306* **** external functions **** 307 real*8 control_fake_mass 308 real*8 ion_amass 309 integer ion_nion 310 external control_fake_mass 311 external ion_amass 312 external ion_nion 313 314 fmass = control_fake_mass() 315 316 317* **** reSet Er0(1) = (1/2)*(g*k*T), where g=number of degrees of freedom **** 318 if (ion_nion().gt.2) then 319 dbl_mb(Er0(1)) = 0.5d0*(3.0d0*dble(ion_nion())-6.0d0)*kb*Tr_new 320 321* **** Dimer molecule. note that the above formula may not work for **** 322* **** linear molecules with more than two atoms. **** 323 else 324 dbl_mb(Er0(1)) = 0.5d0*(1)*kb*Tr_new 325 end if 326 327* **** reSet Er0(2:Nchain) = 1/2*(k*T) **** 328 if (Nchain.gt.1) then 329 call dcopy(Nchain-1,0.5d0*kb*Tr_new,0,dbl_mb(Er0(1)+1),1) 330 end if 331 332* **** total mass **** 333 am = 0.0d0 334 do i=1,ion_nion() 335 am = am + ion_amass(i) 336 end do 337 338c* **** reSet Ee0(1) **** 339c Ee0 = 4.0*kb*Te_new*fmass*dble(ion_nion())/am * eke0_init 340 341* **** Set Ee0(1) - read from rtdb otherwise use current KE **** 342 dbl_mb(Ee0(1))=4.0d0*kb*Te_new*fmass*dble(ion_nion())/am*eke0_init 343 344* **** Set Ee0(2:Mchain) = 1/2*(1/betae), where 1/betae = 2*Ee/Ne **** 345 if (Mchain.gt.1) then 346 betae = dbl_mb(Ee0(1))/Ne_chain 347 call dcopy(Mchain-1,betae,0,dbl_mb(Ee0(1)+1),1) 348 end if 349 350 351 return 352 end 353 354 355 356* *************************** 357* * * 358* * Nose_Newton_Step * 359* * * 360* *************************** 361 362 subroutine Nose_Newton_Step(eke,eki) 363 implicit none 364 real*8 eke,eki 365 366#include "bafdecls.fh" 367#include "nose-hoover.fh" 368 369* **** local variables **** 370 integer m,n 371 real*8 FXe,FXr,dt,a 372 real*8 eke_tmp,ekr_tmp 373 374* **** external functions **** 375 real*8 control_time_step 376 external control_time_step 377 378 dt = control_time_step() 379 380c FXe = 2.0d0*(eke-Ee0) 381c Xe2 = (0.5d0*dt*dt/Qe)*FXe 382 383 eke_tmp = eke 384 do m=1,Mchain-1 385 386* *** integrate thermostat using newton step **** 387 FXe = 2.0d0*(eke_tmp-dbl_mb(Ee0(1)+m-1)) 388 a = dt*(1.0d0 - 0.5d0*dt*dbl_mb(Xem(1)+m)) 389 dbl_mb(Xe2(1)+m-1) = dbl_mb(Xe1(1)+m-1) 390 > + a*dbl_mb(Xem(1)+m-1) 391 > + (0.5d0*dt*dt/dbl_mb(Qe(1)+m-1))*FXe 392 393* **** define kinetic energy for next link in the chain **** 394 eke_tmp = dbl_mb(Xem(1)+m-1) 395 eke_tmp = 0.5d0*dbl_mb(Qe(1)+m-1)*(eke_tmp**2) 396 end do 397 FXe = 2.0d0*(eke_tmp-dbl_mb(Ee0(1)+Mchain-1)) 398 dbl_mb(Xe2(1)+Mchain-1) = dbl_mb(Xe1(1)+Mchain-1) 399 > + dt*dbl_mb(Xem(1)+Mchain-1) 400 > + (0.5d0*dt*dt/dbl_mb(Qe(1)+Mchain-1))*FXe 401 402c FXr = 2.0d0*(eki-Er0) 403c Xr2 = (0.5d0*dt*dt/Qr)*FXr 404 405 ekr_tmp = eki 406 do n=1,Nchain-1 407 408* *** integrate thermostat using newton step **** 409 FXr = 2.0d0*(ekr_tmp-dbl_mb(Er0(1)+n-1)) 410 a = dt*(1.0d0 - 0.5d0*dt*dbl_mb(Xrm(1)+n)) 411 dbl_mb(Xr2(1)+n-1) = dbl_mb(Xr1(1)+n-1) 412 > + a*dbl_mb(Xrm(1)+n-1) 413 > + (0.5d0*dt*dt/dbl_mb(Qr(1)+n-1))*FXr 414 415* **** define kinetic energy for next link in the chain **** 416 ekr_tmp = dbl_mb(Xrm(1)+n-1) 417 ekr_tmp = 0.5d0*dbl_mb(Qr(1)+n-1)*(ekr_tmp**2) 418 end do 419 FXr = 2.0d0*(ekr_tmp-dbl_mb(Er0(1)+Nchain-1)) 420 dbl_mb(Xr2(1)+Nchain-1) = dbl_mb(Xr1(1)+Nchain-1) 421 > + dt*dbl_mb(Xrm(1)+Nchain-1) 422 > + (0.5d0*dt*dt/dbl_mb(Qr(1)+Nchain-1))*FXr 423 424 return 425 end 426 427* *************************** 428* * * 429* * Nose_Verlet_Step * 430* * * 431* *************************** 432 433 subroutine Nose_Verlet_Step(eke,eki) 434 implicit none 435 real*8 eke,eki 436 437#include "bafdecls.fh" 438#include "nose-hoover.fh" 439 440* **** local variables **** 441 integer m,n 442 real*8 eke_tmp,ekr_tmp 443 real*8 FXe,dXe,sse,FXr,dXr,ssr,dt 444 445* **** external functions **** 446 real*8 control_time_step 447 external control_time_step 448 449 dt = control_time_step() 450 451c eke_tmp = eke 452c FXe = 2.0d0*(eke_tmp-Ee0) 453c Xe2 = 2.0d0*Xe1 - Xe0 + (dt*dt/Qe)*FXe 454 455 eke_tmp = eke 456 do m=1,Mchain-1 457* **** define dXe/dt = (3*Xe(t) - 4*Xe(t-dt) + Xe(t-2*dt))/(2*dt) **** 458 dXe = (3.0d0*dbl_mb(Xe1(1)+m) 459 > -4.0d0*dbl_mb(Xe0(1)+m) 460 > + dbl_mb(Xem(1)+m))/(2.0d0*dt) 461 sse = 1.0d0/(1.0d0+0.5d0*dXe*dt) 462 463* *** integrate thermostat using modified verlet **** 464 FXe = 2.0d0*(eke_tmp-dbl_mb(Ee0(1)+m-1)) 465 dbl_mb(Xe2(1)+m-1) = dbl_mb(Xe0(1)+m-1) 466 > + ( dbl_mb(Xe1(1)+m-1) 467 > - dbl_mb(Xe0(1)+m-1) 468 > + (0.5*dt*dt/dbl_mb(Qe(1)+m-1))*FXe 469 > )*2.0d0*sse 470 471* **** define kinetic energy for next link in the chain **** 472 eke_tmp = (dbl_mb(Xe2(1)+m-1)-dbl_mb(Xe0(1)+m-1))/(2.0d0*dt) 473 eke_tmp = 0.5d0*dbl_mb(Qe(1)+m-1)*(eke_tmp**2) 474 end do 475 476* **** Last link of chain **** 477 FXe = 2.0d0*(eke_tmp-dbl_mb(Ee0(1)+Mchain-1)) 478 dbl_mb(Xe2(1)+Mchain-1) = 2.0d0*dbl_mb(Xe1(1)+Mchain-1) 479 > - dbl_mb(Xe0(1)+Mchain-1) 480 > + (dt*dt/dbl_mb(Qe(1)+Mchain-1))*FXe 481 482 483c eki_tmp = eki 484c FXr = 2.0d0*(eki_tmp-Er0) 485c Xr2 = 2.0d0*Xr1 - Xr0 + (dt*dt/Qr)*FXr 486 487 ekr_tmp = eki 488 do n=1,Nchain-1 489* **** define dXe/dt = (3*Xe(t) - 4*Xe(t-dt) + Xe(t-2*dt))/(2*dt) **** 490 dXr = (3.0d0*dbl_mb(Xr1(1)+n) 491 > -4.0d0*dbl_mb(Xr0(1)+n) 492 > + dbl_mb(Xrm(1)+n))/(2.0d0*dt) 493 ssr = 1.0d0/(1.0d0+0.5d0*dXr*dt) 494 495* *** integrate thermostat using modified verlet **** 496 FXr = 2.0d0*(ekr_tmp-dbl_mb(Er0(1)+n-1)) 497 dbl_mb(Xr2(1)+n-1) = dbl_mb(Xr0(1)+n-1) 498 > + ( dbl_mb(Xr1(1)+n-1) 499 > - dbl_mb(Xr0(1)+n-1) 500 > + (0.5*dt*dt/dbl_mb(Qr(1)+n-1))*FXr 501 > )*2.0d0*ssr 502 503* **** define kinetic energy for next link in the chain **** 504 ekr_tmp = (dbl_mb(Xr2(1)+n-1)-dbl_mb(Xr0(1)+n-1))/(2.0d0*dt) 505 ekr_tmp = 0.5d0*dbl_mb(Qr(1)+n-1)*(ekr_tmp**2) 506 end do 507 508* **** Last link of chain **** 509 FXr = 2.0d0*(ekr_tmp-dbl_mb(Er0(1)+Nchain-1)) 510 dbl_mb(Xr2(1)+Nchain-1) = 2.0d0*dbl_mb(Xr1(1)+Nchain-1) 511 > - dbl_mb(Xr0(1)+Nchain-1) 512 > + (dt*dt/dbl_mb(Qr(1)+Nchain-1))*FXr 513 514 515 return 516 end 517 518* *************************** 519* * * 520* * Nose_dXe * 521* * * 522* *************************** 523 524* returns the velocity of the first electronic thermostat 525* used for Newton Step 526 527 real*8 function Nose_dXe() 528 implicit none 529 530#include "bafdecls.fh" 531#include "nose-hoover.fh" 532 533 Nose_dXe = dbl_mb(Xem(1)) 534 return 535 end 536 537* *************************** 538* * * 539* * Nose_dXr * 540* * * 541* *************************** 542 543* returns the velocity of the first ion thermostat 544* used for Newton Step 545 546 real*8 function Nose_dXr() 547 implicit none 548 549#include "bafdecls.fh" 550#include "nose-hoover.fh" 551 552 Nose_dXr = dbl_mb(Xrm(1)) 553 return 554 end 555 556* *************************** 557* * * 558* * Nose_sse * 559* * * 560* *************************** 561 562 real*8 function Nose_sse() 563 implicit none 564 565#include "bafdecls.fh" 566#include "nose-hoover.fh" 567 568 569* ***** local variables **** 570 real*8 dXe,dt 571 572* **** external functions **** 573 real*8 control_time_step 574 external control_time_step 575 576 dt = control_time_step() 577 578 dXe = (3.0d0*dbl_mb(Xe1(1)) 579 > -4.0d0*dbl_mb(Xe0(1)) 580 > + dbl_mb(Xem(1)))/(2.0d0*dt) 581 582 Nose_sse = 1.0d0/(1.0d0+0.5d0*dXe*dt) 583 return 584 end 585 586* *************************** 587* * * 588* * Nose_ssr * 589* * * 590* *************************** 591 592 real*8 function Nose_ssr() 593 implicit none 594 595#include "bafdecls.fh" 596#include "nose-hoover.fh" 597 598* ***** local variables **** 599 real*8 dXr,dt 600 601* **** external functions **** 602 real*8 control_ion_time_step 603 external control_ion_time_step 604 605 dt = control_ion_time_step() 606 607 dXr = (3.0d0*dbl_mb(Xr1(1)) 608 > -4.0d0*dbl_mb(Xr0(1)) 609 > + dbl_mb(Xrm(1)))/(2.0d0*dt) 610 Nose_ssr = 1.0d0/(1.0d0+0.5d0*dXr*dt) 611 return 612 end 613 614* *************************** 615* * * 616* * Nose_shift * 617* * * 618* *************************** 619 620 subroutine Nose_shift() 621 implicit none 622 623#include "bafdecls.fh" 624#include "nose-hoover.fh" 625 626 call dcopy(Mchain,dbl_mb(Xe0(1)),1,dbl_mb(Xem(1)),1) 627 call dcopy(Mchain,dbl_mb(Xe1(1)),1,dbl_mb(Xe0(1)),1) 628 call dcopy(Mchain,dbl_mb(Xe2(1)),1,dbl_mb(Xe1(1)),1) 629 630 call dcopy(Nchain,dbl_mb(Xr0(1)),1,dbl_mb(Xrm(1)),1) 631 call dcopy(Nchain,dbl_mb(Xr1(1)),1,dbl_mb(Xr0(1)),1) 632 call dcopy(Nchain,dbl_mb(Xr2(1)),1,dbl_mb(Xr1(1)),1) 633 return 634 end 635 636* *************************** 637* * * 638* * Nose_e_energy * 639* * * 640* *************************** 641 642 real*8 function Nose_e_energy() 643 implicit none 644 645#include "bafdecls.fh" 646#include "nose-hoover.fh" 647 648* **** local variables **** 649 integer m 650 real*8 dXe,dt,esum 651 652* **** external functions **** 653 real*8 control_time_step 654 external control_time_step 655 656 dt = control_time_step() 657 658 esum = 0.0d0 659 do m=1,Mchain 660 dXe = (3.0d0*dbl_mb(Xe1(1)+m-1) 661 > -4.0d0*dbl_mb(Xe0(1)+m-1) 662 > + dbl_mb(Xem(1)+m-1))/(2.0d0*dt) 663c dXe = (dbl_mb(Xe2(1)+m-1)-dbl_mb(Xe0(1)+m-1))/(2.0d0*dt) 664 665 esum = esum + 0.5d0*dbl_mb(Qe(1)+m-1)*dXe**2 666 > + 2.0d0*dbl_mb(Ee0(1)+m-1)*dbl_mb(Xe1(1)+m-1) 667 end do 668 669 Nose_e_energy = esum 670 return 671 end 672 673* *************************** 674* * * 675* * Nose_r_energy * 676* * * 677* *************************** 678 679 real*8 function Nose_r_energy() 680 implicit none 681 682#include "bafdecls.fh" 683#include "nose-hoover.fh" 684 685* **** local variables **** 686 integer n 687 real*8 dXr,dt,esum 688 689* **** external functions **** 690 real*8 control_time_step 691 external control_time_step 692 693 dt = control_time_step() 694 695 esum = 0.0d0 696 do n=1,Nchain 697 dXr = (3.0d0*dbl_mb(Xr1(1)+n-1) 698 > -4.0d0*dbl_mb(Xr0(1)+n-1) 699 > + dbl_mb(Xrm(1)+n-1))/(2.0d0*dt) 700c dXr = (dbl_mb(Xr2(1)+n-1)-dbl_mb(Xr0(1)+n-1))/(2.0d0*dt) 701 esum = esum + 0.5d0*dbl_mb(Qr(1)+n-1)*dXr**2 702 > + 2.0d0*dbl_mb(Er0(1)+n-1)*dbl_mb(Xr1(1)+n-1) 703 end do 704 705 Nose_r_energy = esum 706 return 707 end 708 709 710* ********************* 711* * * 712* * Nose_Qe * 713* * * 714* ********************* 715 716 real*8 function Nose_Qe(i) 717 implicit none 718 integer i 719 720#include "bafdecls.fh" 721#include "nose-hoover.fh" 722 723 Nose_Qe = dbl_mb(Qe(1)+i-1) 724 return 725 end 726 727 728* ********************* 729* * * 730* * Nose_Pe * 731* * * 732* ********************* 733 734 real*8 function Nose_Pe(i) 735 implicit none 736 integer i 737 738#include "bafdecls.fh" 739#include "nose-hoover.fh" 740 741 Nose_Pe = dbl_mb(Pe(1)+i-1) 742 return 743 end 744 745 746* ********************* 747* * * 748* * Nose_Qr * 749* * * 750* ********************* 751 752 real*8 function Nose_Qr(i) 753 implicit none 754 integer i 755 756#include "bafdecls.fh" 757#include "nose-hoover.fh" 758 759 Nose_Qr = dbl_mb(Qr(1)+i-1) 760 return 761 end 762 763* ********************* 764* * * 765* * Nose_Pr * 766* * * 767* ********************* 768 769 real*8 function Nose_Pr(i) 770 implicit none 771 integer i 772 773#include "bafdecls.fh" 774#include "nose-hoover.fh" 775 776 Nose_Pr = dbl_mb(Pr(1)+i-1) 777 return 778 end 779 780 781* ********************* 782* * * 783* * Nose_Ee0 * 784* * * 785* ********************* 786 787 real*8 function Nose_Ee0(i) 788 implicit none 789 integer i 790 791#include "bafdecls.fh" 792#include "nose-hoover.fh" 793 794 Nose_Ee0 = dbl_mb(Ee0(1)+i-1) 795 return 796 end 797 798* ********************* 799* * * 800* * Nose_Er0 * 801* * * 802* ********************* 803 804 real*8 function Nose_Er0(i) 805 implicit none 806 integer i 807 808#include "bafdecls.fh" 809#include "nose-hoover.fh" 810 811 Nose_Er0 = dbl_mb(Er0(1)+i-1) 812 return 813 end 814 815* ********************* 816* * * 817* * Nose_Mchain * 818* * * 819* ********************* 820 integer function Nose_Mchain() 821 implicit none 822#include "nose-hoover.fh" 823 Nose_Mchain = Mchain 824 return 825 end 826 827* ********************* 828* * * 829* * Nose_Nchain * 830* * * 831* ********************* 832 integer function Nose_Nchain() 833 implicit none 834#include "nose-hoover.fh" 835 Nose_Nchain = Nchain 836 return 837 end 838