TABLE OF CONTENTS
ABINIT/isopress [ Functions ]
NAME
isopress
FUNCTION
performs one half step on isopress parameters according to Martyna et al.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, JCC, JYR, SE) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...) dtion= ionic time step isotemp_data ktemp press= current pressure of the system prtarget= target pressure ucvol= unit cell volume vel= current velocity
OUTPUT
Only updates variables
SIDE EFFECTS
isotemp_data: updates the thermostat parameters (saved variables: bouh !) vel=update the velocities
NOTES
This program is decribed in the following paper Explicit integrators for extended systems dynamics Glenn J Martyna et al. Mol. Phys., 1996, Vol. 87, pp. 1117-1157
PARENTS
pred_isothermal
CHILDREN
dsyev
SOURCE
207 #if defined HAVE_CONFIG_H 208 #include "config.h" 209 #endif 210 211 subroutine isopress(amass,bmass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,qmass,& 212 & strten,strtarget,ucvol,vel,vlogv) 213 214 use m_profiling_abi 215 216 use defs_basis 217 !,isotemp_data) 218 use m_abimover 219 220 !This section has been created automatically by the script Abilint (TD). 221 !Do not modify the following lines by hand. 222 #undef ABI_FUNC 223 #define ABI_FUNC 'isopress' 224 !End of the abilint section 225 226 implicit none 227 228 !Arguments ------------------------------------ 229 !scalars 230 integer,intent(in) :: nnos,natom 231 real(dp),intent(in) :: dtion,ktemp,ucvol,bmass 232 real(dp),intent(inout) :: vlogv 233 real(dp),intent(out) :: ekin 234 type(mttk_type) :: mttk_vars 235 !arrays 236 real(dp),intent(in) :: amass(natom),strtarget(6),strten(6) 237 real(dp),intent(inout) :: vel(3,natom) 238 real(dp),intent(in) :: qmass(:) 239 integer,intent(in) :: iatfix(:,:) 240 241 !Local variables ------------------------------ 242 !scalars 243 integer :: iatom,idir,inos 244 real(dp) :: alocal,glogv,gn1kt,gnkt,nfree,odnf,press,prtarget,scale 245 logical :: DEBUG=.FALSE. 246 !character(len=500) :: message 247 !arrays 248 real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:) 249 250 !*************************************************************************** 251 !Beginning of executable session 252 !*************************************************************************** 253 254 if(DEBUG) then 255 write(std_out,*) ch10,'***',ch10 256 write(std_out,*)' isopress : enter ' 257 write(std_out,*)' strtarget=',strtarget 258 write(std_out,*)' bmass=',bmass 259 write(std_out,*)' dtion=',dtion 260 write(std_out,*)' ktemp=',ktemp 261 write(std_out,*)' natom=',natom 262 write(std_out,*)' nnos=',nnos 263 write(std_out,*)' strten=',strten 264 write(std_out,*)' ucvol=',ucvol 265 end if 266 267 ABI_ALLOCATE(glogs,(nnos)) 268 ABI_ALLOCATE(vlogs,(nnos)) 269 ABI_ALLOCATE(xlogs,(nnos)) 270 glogs(:)=mttk_vars%glogs(:) 271 vlogs(:)=mttk_vars%vlogs(:) 272 xlogs(:)=mttk_vars%xlogs(:) 273 glogv =mttk_vars%glogv 274 scale=one 275 !Compute the ionic kinetic energy 276 nfree=zero 277 ekin=zero 278 do iatom=1,natom 279 do idir=1,3 280 ! Warning : the fixing of atomis is implemented in reduced 281 ! coordinates, so that this expression is wrong 282 if (iatfix(idir,iatom) == 0) then 283 ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2 284 ! Counts the degrees of freedom 285 nfree=nfree+one 286 end if 287 end do 288 end do 289 prtarget=-(strtarget(1)+strtarget(2)+strtarget(3))/three 290 press=-(strten(1)+strten(2)+strten(3))/three 291 gnkt=nfree*ktemp 292 gn1kt=(nfree+one)*ktemp 293 odnf=one+three/nfree 294 !Update the forces 295 glogs(1)=(two*ekin+bmass*vlogv*vlogv-gn1kt)/qmass(1) 296 glogv=(odnf*two*ekin+three*(press-prtarget)*ucvol)/bmass 297 !Update thermostat velocity 298 vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 299 do inos=1,nnos-1 300 alocal=exp(-dtion/eight*vlogs(nnos+1-inos)) 301 vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+& 302 & dtion/four*glogs(nnos-inos)*alocal 303 end do 304 !Update dLog(V)/dt 305 alocal=exp(-dtion/eight*vlogs(1)) 306 vlogv=vlogv*alocal**2+dtion/four*glogv*alocal 307 !Update the particle velocities 308 alocal=exp(-dtion/two*(vlogs(1)+odnf*vlogv)) 309 scale=scale*alocal 310 ekin=ekin*alocal**2 311 glogv=(odnf*two*ekin+three*(press-prtarget)*ucvol)/bmass 312 !Update the thermostat positions 313 do inos=1,nnos 314 xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two 315 end do 316 !Update dLog(V)/dt 317 alocal=exp(-dtion/eight*vlogs(1)) 318 vlogv=vlogv*alocal**2+dtion/four*glogv*alocal 319 !Update the forces 320 glogs(1)=(two*ekin+bmass*vlogv*vlogv-gn1kt)/qmass(1) 321 !Update the thermostat velocities 322 do inos=1,nnos-1 323 alocal=exp(-dtion/eight*vlogs(inos+1)) 324 vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal 325 glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1) 326 end do 327 vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 328 vel(:,:)=vel(:,:)*scale 329 !Compute the ionic kinetic energy 330 ekin=zero 331 do iatom=1,natom 332 do idir=1,3 333 ! Warning : the fixing of atomis is implemented in reduced 334 ! coordinates, so that this expression is wrong 335 if (iatfix(idir,iatom) == 0) then 336 ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2 337 end if 338 end do 339 end do 340 !Compute the thermostat kinetic energy and add it to the ionic one 341 !First thermostat 342 ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*(nfree+one)*ktemp 343 !Other thermostats 344 do inos=2,nnos 345 ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp 346 end do 347 !Barostat 348 ekin=ekin+half*bmass*vlogv**2+prtarget*ucvol 349 ABI_DEALLOCATE(glogs) 350 ABI_DEALLOCATE(vlogs) 351 ABI_DEALLOCATE(xlogs) 352 353 !DEBUG 354 !write(std_out,*) 'EKIN',ekin 355 !write(std_out,*) 'VLOGV',vlogv 356 !write(std_out,*)'ekin added T',half*qmass(:)*vlogs(:)**2,xlogs(:)*(nfree)*ktemp 357 !write(std_out,*)'ekin added P',half*bmass*vlogv**2,prtarget*ucvol 358 !write(std_out,*)'ekin last',ekin 359 !ENDDEBUG 360 end subroutine isopress
ABINIT/isostress [ Functions ]
NAME
isostress
FUNCTION
performs one half step on isostress parameters according to Martyna et al.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, JCC, JYR, SE) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...) dtion= ionic time step isotemp_data ktemp press= current pressure of the system prtarget= target pressure ucvol= unit cell volume vel= current velocity
OUTPUT
Only updates variables
SIDE EFFECTS
isotemp_data: updates the thermostat parameters (saved variables: bouh !) vel=update the velocities
NOTES
This program is decribed in the following paper Explicit integrators for extended systems dynamics Glenn J Martyna et al. Mol. Phys., 1996, Vol. 87, pp. 1117-1157
PARENTS
pred_isothermal
CHILDREN
dsyev
SOURCE
405 #if defined HAVE_CONFIG_H 406 #include "config.h" 407 #endif 408 409 subroutine isostress(amass,bmass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,& 410 & qmass,strten,strtarget,ucvol,vel) 411 412 use m_profiling_abi 413 414 use defs_basis 415 !,isotemp_data) 416 use m_abimover 417 use m_linalg_interfaces 418 419 !This section has been created automatically by the script Abilint (TD). 420 !Do not modify the following lines by hand. 421 #undef ABI_FUNC 422 #define ABI_FUNC 'isostress' 423 !End of the abilint section 424 425 implicit none 426 427 !Arguments ------------------------------------ 428 !scalars 429 integer,intent(in) :: natom,nnos 430 real(dp),intent(in) :: dtion,ktemp,ucvol,bmass 431 real(dp),intent(out) :: ekin 432 type(mttk_type) :: mttk_vars 433 !arrays 434 real(dp),intent(in) :: amass(natom),strtarget(6),strten(6),qmass(:) 435 real(dp),intent(inout) :: vel(3,natom) 436 integer, intent(in) :: iatfix(:,:) 437 438 !Local variables ------------------------------ 439 !scalars 440 integer,parameter :: lwork=8 441 integer :: iatom,idir,info,inos,jdir 442 real(dp) :: akinb,alocal,gn1kt,gnd2kt,nfree,odnf,trvg !,gnkt,scale 443 logical :: DEBUG=.FALSE. 444 !character(len=500) :: message 445 !arrays 446 real(dp) :: akin(3,3),expdiag(3),gboxg(3,3),identity(3,3),press(3,3) 447 real(dp) :: prtarget(3,3),tvtemp(3,3),uv(3),vboxg(3,3),veig(3),vtemp(3,3) 448 real(dp) :: work(lwork) 449 real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:) 450 451 !*************************************************************************** 452 !Beginning of executable session 453 !*************************************************************************** 454 455 if(DEBUG) then 456 write(std_out,*) ch10,'***',ch10 457 write(std_out,*)' isostress : enter ' 458 write(std_out,*)' strtarget=',strtarget 459 write(std_out,*)' bmass=',bmass 460 write(std_out,*)' dtion=',dtion 461 write(std_out,*)' ktemp=',ktemp 462 write(std_out,*)' natom=',natom 463 write(std_out,*)' nnos=',nnos 464 write(std_out,*)' strten=',strten 465 write(std_out,*)' ucvol=',ucvol 466 end if 467 468 ABI_ALLOCATE(glogs,(nnos)) 469 ABI_ALLOCATE(vlogs,(nnos)) 470 ABI_ALLOCATE(xlogs,(nnos)) 471 glogs(:)=mttk_vars%glogs(:) 472 vlogs(:)=mttk_vars%vlogs(:) 473 xlogs(:)=mttk_vars%xlogs(:) 474 vboxg(:,:)=mttk_vars%vboxg(:,:) 475 identity(:,:)=zero 476 do idir=1,3 477 identity(idir,idir)=one 478 end do 479 480 !write(std_out,*) 'isostress 02' 481 !########################################################## 482 !### 03. Compute the ionic kinetic energy 483 484 nfree=zero 485 ekin=zero 486 do iatom=1,natom 487 do idir=1,3 488 ! Warning : the fixing of atomis is implemented in reduced 489 ! coordinates, so that this exprtargetion is wrong 490 if (iatfix(idir,iatom) == 0) then 491 ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2 492 ! Counts the degrees of freedom 493 nfree=nfree+one 494 end if 495 end do 496 end do 497 498 gn1kt=(nfree+one)*ktemp 499 gnd2kt=(nfree+9)*ktemp 500 odnf=one+three/nfree 501 akin(:,:)=zero 502 do iatom=1,natom 503 do idir=1,3 504 do jdir=1,3 505 ! Warning : the fixing of atomis is implemented in reduced 506 ! coordinates, so that this expression is wrong 507 akin(idir,jdir)=akin(idir,jdir)+0.5d0*amass(iatom)*vel(idir,iatom)*vel(jdir,iatom) 508 end do 509 end do 510 end do 511 akinb=zero 512 do idir=1,3 513 do jdir=1,3 514 akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2 515 end do 516 end do 517 !Compute the pressure: from Voigt to tensor notation+kinetic energy 518 do idir=1,3 519 press(idir,idir)=-strten(idir) 520 prtarget(idir,idir)=-strtarget(idir) 521 end do 522 press(3,2)=-strten(4); press(1,3)=-strten(5); press(2,1)=-strten(6) 523 prtarget(3,2)=-strtarget(4); prtarget(1,3)=-strtarget(5); prtarget(2,1)=-strtarget(6) 524 press(2,3)=press(3,2); press(3,1)=press(1,3); press(1,2)=press(2,1) 525 prtarget(2,3)=prtarget(3,2); prtarget(3,1)=prtarget(1,3); prtarget(1,2)=prtarget(2,1) 526 !Update the forces 527 glogs(1)=(two*ekin+two*akinb-gnd2kt)/qmass(1) 528 gboxg(:,:)=(two*ekin/nfree*identity(:,:)+two*akin(:,:)+(press(:,:)-prtarget(:,:))*ucvol)/bmass 529 !Update thermostat velocity 530 if (nnos > 0) vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 531 do inos=1,nnos-1 532 alocal=exp(-dtion/eight*vlogs(nnos+1-inos)) 533 vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+& 534 & dtion/four*glogs(nnos-inos)*alocal 535 end do 536 !Update box velocity 537 alocal=exp(-dtion/eight*vlogs(1)) 538 539 if(DEBUG) then 540 write(std_out,*)' gboxg(:,:)=',gboxg(:,:) 541 write(std_out,*)' vboxg(:,:)=',vboxg(:,:) 542 write(std_out,*)' alocal=',alocal 543 end if 544 545 vboxg(:,:)=vboxg(:,:)*alocal**2+dtion/four*gboxg(:,:)*alocal 546 !Update the thermostat positions 547 do inos=1,nnos 548 xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two 549 end do 550 !Update the particle velocities 551 trvg=(vboxg(1,1)+vboxg(2,2)+vboxg(3,3))/nfree 552 vtemp(:,:)=vboxg(:,:)+(trvg+vlogs(1))*identity(:,:) 553 call dsyev('V','U',3,vtemp,3,veig,work,lwork,info) 554 !On exit, we have vtemp=U such that tU vtemp U = veig 555 tvtemp(:,:)=transpose(vtemp) 556 557 if(DEBUG) then 558 write(std_out,*)' vboxg(:,:)=',vboxg(:,:) 559 write(std_out,*)' vtemp(:,:)=',vtemp(:,:) 560 write(std_out,*)' veig(:)=',veig(:) 561 end if 562 563 expdiag(1)=exp(-veig(1)*dtion/two) 564 expdiag(2)=exp(-veig(2)*dtion/two) 565 expdiag(3)=exp(-veig(3)*dtion/two) 566 if(DEBUG) then 567 write(std_out,*)' isostress : expdiag(:)=',expdiag(:) ! Do not remove this line : seems to be needed for g95 compilo 568 end if 569 do iatom=1,natom 570 uv(:)=matmul(tvtemp,vel(:,iatom)) 571 uv(:)=uv(:)*expdiag(:) 572 vel(:,iatom)=matmul(vtemp,uv) 573 end do 574 !Compute the ionic kinetic energy 575 nfree=zero 576 ekin=zero 577 do iatom=1,natom 578 do idir=1,3 579 ! Warning : the fixing of atomis is implemented in reduced 580 ! coordinates, so that this expression is wrong 581 if (iatfix(idir,iatom) == 0) then 582 ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2 583 if(DEBUG) then 584 write(std_out,*)'kin',iatom,ekin,vel(idir,iatom) 585 end if 586 ! Counts the degrees of freedom 587 nfree=nfree+one 588 end if 589 end do 590 end do 591 gn1kt=(nfree+one)*ktemp 592 gnd2kt=(nfree+9)*ktemp 593 odnf=one+three/nfree 594 akin(:,:)=zero 595 do iatom=1,natom 596 do idir=1,3 597 do jdir=1,3 598 ! Warning : the fixing of atomis is implemented in reduced 599 ! coordinates, so that this expression is wrong 600 akin(idir,jdir)=akin(idir,jdir)+0.5d0*amass(iatom)*vel(idir,iatom)*vel(jdir,iatom) 601 end do 602 end do 603 end do 604 gboxg(:,:)=(two*ekin/nfree*identity(:,:)+two*akin(:,:)+(press(:,:)-prtarget(:,:))*ucvol)/bmass 605 !Update box velocity 606 alocal=exp(-dtion/eight*vlogs(1)) 607 vboxg(:,:)=vboxg(:,:)*alocal**2+dtion/four*gboxg(:,:)*alocal 608 !Compute the box kinetic energy 609 akinb=zero 610 do idir=1,3 611 do jdir=1,3 612 akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2 613 end do 614 end do 615 glogs(1)=(two*ekin+two*akinb-gnd2kt)/qmass(1) 616 !Update the thermostat velocities 617 do inos=1,nnos-1 618 alocal=exp(-dtion/eight*vlogs(inos+1)) 619 vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal 620 glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1) 621 end do 622 if (nnos > 0) vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 623 !Compute the ionic kinetic energy 624 ekin=zero 625 do iatom=1,natom 626 do idir=1,3 627 ! Warning : the fixing of atomis is implemented in reduced 628 ! coordinates, so that this expression is wrong 629 if (iatfix(idir,iatom) == 0) then 630 ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2 631 end if 632 end do 633 end do 634 !Compute the thermostat kinetic energy and add it to the ionic one 635 !First thermostat 636 ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*(nfree+nine)*ktemp 637 !Other thermostats 638 do inos=2,nnos 639 ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp 640 end do 641 !Barostat kinetic energy 642 akinb=zero 643 do idir=1,3 644 do jdir=1,3 645 akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2 646 end do 647 end do 648 !ekin is the invariant minus the potential energy 649 ekin=ekin+akinb+prtarget(1,1)*ucvol 650 651 mttk_vars%vboxg(:,:)=vboxg(:,:) 652 mttk_vars%glogs(:)=glogs(:) 653 mttk_vars%vlogs(:)=vlogs(:) 654 mttk_vars%xlogs(:)=xlogs(:) 655 ABI_DEALLOCATE(glogs) 656 ABI_DEALLOCATE(vlogs) 657 ABI_DEALLOCATE(xlogs) 658 659 if(DEBUG) then 660 ! write(std_out,*)'ekin added T',half*qmass(:)*vlogs(:)**2,xlogs(:)*(nfree)*ktemp 661 ! write(std_out,*)'ekin added P',akinb,prtarget*ucvol 662 ! write(std_out,*)'ekin last',ekin 663 write(std_out,*) ch10,' exiting from isostress',ch10 664 end if 665 666 end subroutine isostress
ABINIT/isotemp [ Functions ]
NAME
isotemp
FUNCTION
performs one half step on isotemp parameters according to Martyna et al.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, JCC, JYR, SE) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...) dtion= isotemp_data ktemp vel
OUTPUT
Only updates variables
SIDE EFFECTS
isotemp_data: updates the thermostat parameters vel=update the velocities
NOTES
This program is decribed in the following paper Explicit integrators for extended systems dynamics Glenn J Martyna et al. Mol. Phys., 1996, Vol. 87, pp. 1117-1157
PARENTS
pred_isothermal
CHILDREN
dsyev
SOURCE
43 #if defined HAVE_CONFIG_H 44 #include "config.h" 45 #endif 46 47 #include "abi_common.h" 48 49 subroutine isotemp(amass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,qmass,vel) 50 51 use m_profiling_abi 52 53 use defs_basis 54 use m_abimover 55 56 !This section has been created automatically by the script Abilint (TD). 57 !Do not modify the following lines by hand. 58 #undef ABI_FUNC 59 #define ABI_FUNC 'isotemp' 60 !End of the abilint section 61 62 implicit none 63 64 !Arguments ------------------------------------ 65 !scalars 66 integer,intent(in) :: natom 67 integer,intent(in) :: nnos 68 real(dp),intent(in) :: dtion,ktemp 69 real(dp),intent(out) :: ekin 70 type(mttk_type) :: mttk_vars 71 !arrays 72 real(dp),intent(in) :: amass(natom) 73 real(dp),intent(inout) :: vel(3,natom) 74 real(dp),intent(in) :: qmass(:) 75 integer,intent(in) :: iatfix(:,:) 76 77 !Local variables ------------------------------ 78 !scalars 79 integer :: iatom,idir,inos 80 real(dp) :: alocal,gnkt,nfree,scale 81 !character(len=500) :: message 82 !arrays 83 real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:) 84 85 !*************************************************************************** 86 !Beginning of executable session 87 !*************************************************************************** 88 89 ABI_ALLOCATE(glogs,(nnos)) 90 ABI_ALLOCATE(vlogs,(nnos)) 91 ABI_ALLOCATE(xlogs,(nnos)) 92 glogs(:)=mttk_vars%glogs(:) 93 vlogs(:)=mttk_vars%vlogs(:) 94 xlogs(:)=mttk_vars%xlogs(:) 95 scale=one 96 !Compute the ionic kinetic energy 97 nfree=zero 98 ekin=zero 99 do iatom=1,natom 100 do idir=1,3 101 ! Warning : the fixing of atomis is implemented in reduced 102 ! coordinates, so that this expression is wrong 103 if (iatfix(idir,iatom) == 0) then 104 ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2 105 ! Counts the degrees of freedom 106 nfree=nfree+one 107 end if 108 end do 109 end do 110 gnkt=nfree*ktemp 111 !Update the forces 112 glogs(1)=(two*ekin-gnkt)/qmass(1) 113 vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 114 do inos=1,nnos-1 115 alocal=exp(-dtion/eight*vlogs(nnos+1-inos)) 116 vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+& 117 & dtion/four*glogs(nnos-inos)*alocal 118 end do 119 !Update the particle velocities 120 alocal=exp(-dtion/two*vlogs(1)) 121 scale=scale*alocal 122 !Update the forces 123 glogs(1)=(scale*scale*two*ekin-gnkt)/qmass(1) 124 !Update the thermostat positions 125 do inos=1,nnos 126 xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two 127 end do 128 !Update the thermostat velocities 129 do inos=1,nnos-1 130 alocal=exp(-dtion/eight*vlogs(inos+1)) 131 vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal 132 glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1) 133 end do 134 vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 135 vel(:,:)=vel(:,:)*scale 136 !Compute the ionic kinetic energy 137 ekin=zero 138 do iatom=1,natom 139 do idir=1,3 140 ! Warning : the fixing of atomis is implemented in reduced 141 ! coordinates, so that this expression is wrong 142 if (iatfix(idir,iatom) == 0) then 143 ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2 144 end if 145 end do 146 end do 147 !Compute the thermostat kinetic energy and add it to the ionic one 148 ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*nfree*ktemp 149 do inos=2,nnos 150 ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp 151 end do 152 mttk_vars%glogs(:)=glogs(:) 153 mttk_vars%vlogs(:)=vlogs(:) 154 mttk_vars%xlogs(:)=xlogs(:) 155 ABI_DEALLOCATE(glogs) 156 ABI_DEALLOCATE(vlogs) 157 ABI_DEALLOCATE(xlogs) 158 !DEBUG 159 !write(std_out,*)'ekin added',half*qmass(1)*vlogs(1)**2,xlogs(1)*(nfree)*ktemp 160 !ENDEBUG 161 162 end subroutine isotemp