TABLE OF CONTENTS
ABINIT/isopress [ Functions ]
NAME
isopress
FUNCTION
performs one half step on isopress parameters according to Martyna et al.
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
SOURCE
746 subroutine isopress(amass,bmass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,qmass,& 747 & strten,strtarget,ucvol,vel,vlogv) 748 749 implicit none 750 751 !Arguments ------------------------------------ 752 !scalars 753 integer,intent(in) :: nnos,natom 754 real(dp),intent(in) :: dtion,ktemp,ucvol,bmass 755 real(dp),intent(inout) :: vlogv 756 real(dp),intent(out) :: ekin 757 type(mttk_type) :: mttk_vars 758 !arrays 759 real(dp),intent(in) :: amass(natom),strtarget(6),strten(6) 760 real(dp),intent(inout) :: vel(3,natom) 761 real(dp),intent(in) :: qmass(:) 762 integer,intent(in) :: iatfix(:,:) 763 764 !Local variables ------------------------------ 765 !scalars 766 integer :: iatom,idir,inos 767 real(dp) :: alocal,glogv,gn1kt,gnkt,nfree,odnf,press,prtarget,scale 768 logical :: DEBUG=.FALSE. 769 !character(len=500) :: message 770 !arrays 771 real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:) 772 773 !*************************************************************************** 774 !Beginning of executable session 775 !*************************************************************************** 776 777 if(DEBUG) then 778 write(std_out,*) ch10,'***',ch10 779 write(std_out,*)' isopress : enter ' 780 write(std_out,*)' strtarget=',strtarget 781 write(std_out,*)' bmass=',bmass 782 write(std_out,*)' dtion=',dtion 783 write(std_out,*)' ktemp=',ktemp 784 write(std_out,*)' natom=',natom 785 write(std_out,*)' nnos=',nnos 786 write(std_out,*)' strten=',strten 787 write(std_out,*)' ucvol=',ucvol 788 end if 789 790 ABI_MALLOC(glogs,(nnos)) 791 ABI_MALLOC(vlogs,(nnos)) 792 ABI_MALLOC(xlogs,(nnos)) 793 glogs(:)=mttk_vars%glogs(:) 794 vlogs(:)=mttk_vars%vlogs(:) 795 xlogs(:)=mttk_vars%xlogs(:) 796 glogv =mttk_vars%glogv 797 scale=one 798 !Compute the ionic kinetic energy 799 nfree=zero 800 ekin=zero 801 do iatom=1,natom 802 do idir=1,3 803 ! Warning : the fixing of atomis is implemented in reduced 804 ! coordinates, so that this expression is wrong 805 if (iatfix(idir,iatom) == 0) then 806 ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2 807 ! Counts the degrees of freedom 808 nfree=nfree+one 809 end if 810 end do 811 end do 812 prtarget=-(strtarget(1)+strtarget(2)+strtarget(3))/three 813 press=-(strten(1)+strten(2)+strten(3))/three 814 gnkt=nfree*ktemp 815 gn1kt=(nfree+one)*ktemp 816 odnf=one+three/nfree 817 !Update the forces 818 glogs(1)=(two*ekin+bmass*vlogv*vlogv-gn1kt)/qmass(1) 819 glogv=(odnf*two*ekin+three*(press-prtarget)*ucvol)/bmass 820 !Update thermostat velocity 821 vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 822 do inos=1,nnos-1 823 alocal=exp(-dtion/eight*vlogs(nnos+1-inos)) 824 vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+& 825 & dtion/four*glogs(nnos-inos)*alocal 826 end do 827 !Update dLog(V)/dt 828 alocal=exp(-dtion/eight*vlogs(1)) 829 vlogv=vlogv*alocal**2+dtion/four*glogv*alocal 830 !Update the particle velocities 831 alocal=exp(-dtion/two*(vlogs(1)+odnf*vlogv)) 832 scale=scale*alocal 833 ekin=ekin*alocal**2 834 glogv=(odnf*two*ekin+three*(press-prtarget)*ucvol)/bmass 835 !Update the thermostat positions 836 do inos=1,nnos 837 xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two 838 end do 839 !Update dLog(V)/dt 840 alocal=exp(-dtion/eight*vlogs(1)) 841 vlogv=vlogv*alocal**2+dtion/four*glogv*alocal 842 !Update the forces 843 glogs(1)=(two*ekin+bmass*vlogv*vlogv-gn1kt)/qmass(1) 844 !Update the thermostat velocities 845 do inos=1,nnos-1 846 alocal=exp(-dtion/eight*vlogs(inos+1)) 847 vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal 848 glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1) 849 end do 850 vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 851 vel(:,:)=vel(:,:)*scale 852 !Compute the ionic kinetic energy 853 ekin=zero 854 do iatom=1,natom 855 do idir=1,3 856 ! Warning : the fixing of atomis is implemented in reduced 857 ! coordinates, so that this expression is wrong 858 if (iatfix(idir,iatom) == 0) then 859 ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2 860 end if 861 end do 862 end do 863 !Compute the thermostat kinetic energy and add it to the ionic one 864 !First thermostat 865 ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*(nfree+one)*ktemp 866 !Other thermostats 867 do inos=2,nnos 868 ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp 869 end do 870 !Barostat 871 ekin=ekin+half*bmass*vlogv**2+prtarget*ucvol 872 ABI_FREE(glogs) 873 ABI_FREE(vlogs) 874 ABI_FREE(xlogs) 875 876 !DEBUG 877 !write(std_out,*) 'EKIN',ekin 878 !write(std_out,*) 'VLOGV',vlogv 879 !write(std_out,*)'ekin added T',half*qmass(:)*vlogs(:)**2,xlogs(:)*(nfree)*ktemp 880 !write(std_out,*)'ekin added P',half*bmass*vlogv**2,prtarget*ucvol 881 !write(std_out,*)'ekin last',ekin 882 !ENDDEBUG 883 end subroutine isopress
ABINIT/isostress [ Functions ]
NAME
isostress
FUNCTION
performs one half step on isostress parameters according to Martyna et al.
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
SOURCE
914 subroutine isostress(amass,bmass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,& 915 & qmass,strten,strtarget,ucvol,vel) 916 917 use m_linalg_interfaces 918 implicit none 919 920 !Arguments ------------------------------------ 921 !scalars 922 integer,intent(in) :: natom,nnos 923 real(dp),intent(in) :: dtion,ktemp,ucvol,bmass 924 real(dp),intent(out) :: ekin 925 type(mttk_type) :: mttk_vars 926 !arrays 927 real(dp),intent(in) :: amass(natom),strtarget(6),strten(6),qmass(:) 928 real(dp),intent(inout) :: vel(3,natom) 929 integer, intent(in) :: iatfix(:,:) 930 931 !Local variables ------------------------------ 932 !scalars 933 integer,parameter :: lwork=8 934 integer :: iatom,idir,info,inos,jdir 935 real(dp) :: akinb,alocal,gn1kt,gnd2kt,nfree,odnf,trvg !,gnkt,scale 936 logical :: DEBUG=.FALSE. 937 !character(len=500) :: message 938 !arrays 939 real(dp) :: akin(3,3),expdiag(3),gboxg(3,3),identity(3,3),press(3,3) 940 real(dp) :: prtarget(3,3),tvtemp(3,3),uv(3),vboxg(3,3),veig(3),vtemp(3,3) 941 real(dp) :: work(lwork) 942 real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:) 943 944 !*************************************************************************** 945 !Beginning of executable session 946 !*************************************************************************** 947 948 if(DEBUG) then 949 write(std_out,*) ch10,'***',ch10 950 write(std_out,*)' isostress : enter ' 951 write(std_out,*)' strtarget=',strtarget 952 write(std_out,*)' bmass=',bmass 953 write(std_out,*)' dtion=',dtion 954 write(std_out,*)' ktemp=',ktemp 955 write(std_out,*)' natom=',natom 956 write(std_out,*)' nnos=',nnos 957 write(std_out,*)' strten=',strten 958 write(std_out,*)' ucvol=',ucvol 959 end if 960 961 ABI_MALLOC(glogs,(nnos)) 962 ABI_MALLOC(vlogs,(nnos)) 963 ABI_MALLOC(xlogs,(nnos)) 964 glogs(:)=mttk_vars%glogs(:) 965 vlogs(:)=mttk_vars%vlogs(:) 966 xlogs(:)=mttk_vars%xlogs(:) 967 vboxg(:,:)=mttk_vars%vboxg(:,:) 968 identity(:,:)=zero 969 do idir=1,3 970 identity(idir,idir)=one 971 end do 972 973 !write(std_out,*) 'isostress 02' 974 !########################################################## 975 !### 03. Compute the ionic kinetic energy 976 977 nfree=zero 978 ekin=zero 979 do iatom=1,natom 980 do idir=1,3 981 ! Warning : the fixing of atomis is implemented in reduced 982 ! coordinates, so that this exprtargetion is wrong 983 if (iatfix(idir,iatom) == 0) then 984 ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2 985 ! Counts the degrees of freedom 986 nfree=nfree+one 987 end if 988 end do 989 end do 990 991 gn1kt=(nfree+one)*ktemp 992 gnd2kt=(nfree+9)*ktemp 993 odnf=one+three/nfree 994 akin(:,:)=zero 995 do iatom=1,natom 996 do idir=1,3 997 do jdir=1,3 998 ! Warning : the fixing of atomis is implemented in reduced 999 ! coordinates, so that this expression is wrong 1000 akin(idir,jdir)=akin(idir,jdir)+0.5d0*amass(iatom)*vel(idir,iatom)*vel(jdir,iatom) 1001 end do 1002 end do 1003 end do 1004 akinb=zero 1005 do idir=1,3 1006 do jdir=1,3 1007 akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2 1008 end do 1009 end do 1010 !Compute the pressure: from Voigt to tensor notation+kinetic energy 1011 do idir=1,3 1012 press(idir,idir)=-strten(idir) 1013 prtarget(idir,idir)=-strtarget(idir) 1014 end do 1015 press(3,2)=-strten(4); press(1,3)=-strten(5); press(2,1)=-strten(6) 1016 prtarget(3,2)=-strtarget(4); prtarget(1,3)=-strtarget(5); prtarget(2,1)=-strtarget(6) 1017 press(2,3)=press(3,2); press(3,1)=press(1,3); press(1,2)=press(2,1) 1018 prtarget(2,3)=prtarget(3,2); prtarget(3,1)=prtarget(1,3); prtarget(1,2)=prtarget(2,1) 1019 !Update the forces 1020 glogs(1)=(two*ekin+two*akinb-gnd2kt)/qmass(1) 1021 gboxg(:,:)=(two*ekin/nfree*identity(:,:)+two*akin(:,:)+(press(:,:)-prtarget(:,:))*ucvol)/bmass 1022 !Update thermostat velocity 1023 if (nnos > 0) vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 1024 do inos=1,nnos-1 1025 alocal=exp(-dtion/eight*vlogs(nnos+1-inos)) 1026 vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+& 1027 & dtion/four*glogs(nnos-inos)*alocal 1028 end do 1029 !Update box velocity 1030 alocal=exp(-dtion/eight*vlogs(1)) 1031 1032 if(DEBUG) then 1033 write(std_out,*)' gboxg(:,:)=',gboxg(:,:) 1034 write(std_out,*)' vboxg(:,:)=',vboxg(:,:) 1035 write(std_out,*)' alocal=',alocal 1036 end if 1037 1038 vboxg(:,:)=vboxg(:,:)*alocal**2+dtion/four*gboxg(:,:)*alocal 1039 !Update the thermostat positions 1040 do inos=1,nnos 1041 xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two 1042 end do 1043 !Update the particle velocities 1044 trvg=(vboxg(1,1)+vboxg(2,2)+vboxg(3,3))/nfree 1045 vtemp(:,:)=vboxg(:,:)+(trvg+vlogs(1))*identity(:,:) 1046 call dsyev('V','U',3,vtemp,3,veig,work,lwork,info) 1047 !On exit, we have vtemp=U such that tU vtemp U = veig 1048 tvtemp(:,:)=transpose(vtemp) 1049 1050 if(DEBUG) then 1051 write(std_out,*)' vboxg(:,:)=',vboxg(:,:) 1052 write(std_out,*)' vtemp(:,:)=',vtemp(:,:) 1053 write(std_out,*)' veig(:)=',veig(:) 1054 end if 1055 1056 expdiag(1)=exp(-veig(1)*dtion/two) 1057 expdiag(2)=exp(-veig(2)*dtion/two) 1058 expdiag(3)=exp(-veig(3)*dtion/two) 1059 !if(DEBUG) then 1060 ! write(std_out,*)' isostress : expdiag(:)=',expdiag(:) ! Do not remove this line : seems to be needed for g95 compilo 1061 !end if 1062 do iatom=1,natom 1063 uv(:)=matmul(tvtemp,vel(:,iatom)) 1064 uv(:)=uv(:)*expdiag(:) 1065 vel(:,iatom)=matmul(vtemp,uv) 1066 end do 1067 !Compute the ionic kinetic energy 1068 nfree=zero 1069 ekin=zero 1070 do iatom=1,natom 1071 do idir=1,3 1072 ! Warning : the fixing of atomis is implemented in reduced 1073 ! coordinates, so that this expression is wrong 1074 if (iatfix(idir,iatom) == 0) then 1075 ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2 1076 if(DEBUG) then 1077 write(std_out,*)'kin',iatom,ekin,vel(idir,iatom) 1078 end if 1079 ! Counts the degrees of freedom 1080 nfree=nfree+one 1081 end if 1082 end do 1083 end do 1084 gn1kt=(nfree+one)*ktemp 1085 gnd2kt=(nfree+9)*ktemp 1086 odnf=one+three/nfree 1087 akin(:,:)=zero 1088 do iatom=1,natom 1089 do idir=1,3 1090 do jdir=1,3 1091 ! Warning : the fixing of atomis is implemented in reduced 1092 ! coordinates, so that this expression is wrong 1093 akin(idir,jdir)=akin(idir,jdir)+0.5d0*amass(iatom)*vel(idir,iatom)*vel(jdir,iatom) 1094 end do 1095 end do 1096 end do 1097 gboxg(:,:)=(two*ekin/nfree*identity(:,:)+two*akin(:,:)+(press(:,:)-prtarget(:,:))*ucvol)/bmass 1098 !Update box velocity 1099 alocal=exp(-dtion/eight*vlogs(1)) 1100 vboxg(:,:)=vboxg(:,:)*alocal**2+dtion/four*gboxg(:,:)*alocal 1101 !Compute the box kinetic energy 1102 akinb=zero 1103 do idir=1,3 1104 do jdir=1,3 1105 akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2 1106 end do 1107 end do 1108 glogs(1)=(two*ekin+two*akinb-gnd2kt)/qmass(1) 1109 !Update the thermostat velocities 1110 do inos=1,nnos-1 1111 alocal=exp(-dtion/eight*vlogs(inos+1)) 1112 vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal 1113 glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1) 1114 end do 1115 if (nnos > 0) vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 1116 !Compute the ionic kinetic energy 1117 ekin=zero 1118 do iatom=1,natom 1119 do idir=1,3 1120 ! Warning : the fixing of atomis is implemented in reduced 1121 ! coordinates, so that this expression is wrong 1122 if (iatfix(idir,iatom) == 0) then 1123 ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2 1124 end if 1125 end do 1126 end do 1127 !Compute the thermostat kinetic energy and add it to the ionic one 1128 !First thermostat 1129 ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*(nfree+nine)*ktemp 1130 !Other thermostats 1131 do inos=2,nnos 1132 ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp 1133 end do 1134 !Barostat kinetic energy 1135 akinb=zero 1136 do idir=1,3 1137 do jdir=1,3 1138 akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2 1139 end do 1140 end do 1141 !ekin is the invariant minus the potential energy 1142 ekin=ekin+akinb+prtarget(1,1)*ucvol 1143 1144 mttk_vars%vboxg(:,:)=vboxg(:,:) 1145 mttk_vars%glogs(:)=glogs(:) 1146 mttk_vars%vlogs(:)=vlogs(:) 1147 mttk_vars%xlogs(:)=xlogs(:) 1148 ABI_FREE(glogs) 1149 ABI_FREE(vlogs) 1150 ABI_FREE(xlogs) 1151 1152 if(DEBUG) then 1153 ! write(std_out,*)'ekin added T',half*qmass(:)*vlogs(:)**2,xlogs(:)*(nfree)*ktemp 1154 ! write(std_out,*)'ekin added P',akinb,prtarget*ucvol 1155 ! write(std_out,*)'ekin last',ekin 1156 write(std_out,*) ch10,' exiting from isostress',ch10 1157 end if 1158 1159 end subroutine isostress
ABINIT/isotemp [ Functions ]
NAME
isotemp
FUNCTION
performs one half step on isotemp parameters according to Martyna et al.
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
SOURCE
613 subroutine isotemp(amass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,qmass,vel) 614 615 implicit none 616 617 !Arguments ------------------------------------ 618 !scalars 619 integer,intent(in) :: natom 620 integer,intent(in) :: nnos 621 real(dp),intent(in) :: dtion,ktemp 622 real(dp),intent(out) :: ekin 623 type(mttk_type) :: mttk_vars 624 !arrays 625 real(dp),intent(in) :: amass(natom) 626 real(dp),intent(inout) :: vel(3,natom) 627 real(dp),intent(in) :: qmass(:) 628 integer,intent(in) :: iatfix(:,:) 629 630 !Local variables ------------------------------ 631 !scalars 632 integer :: iatom,idir,inos 633 real(dp) :: alocal,gnkt,nfree,scale 634 !character(len=500) :: message 635 !arrays 636 real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:) 637 638 !*************************************************************************** 639 !Beginning of executable session 640 !*************************************************************************** 641 642 ABI_MALLOC(glogs,(nnos)) 643 ABI_MALLOC(vlogs,(nnos)) 644 ABI_MALLOC(xlogs,(nnos)) 645 glogs(:)=mttk_vars%glogs(:) 646 vlogs(:)=mttk_vars%vlogs(:) 647 xlogs(:)=mttk_vars%xlogs(:) 648 scale=one 649 !Compute the ionic kinetic energy 650 nfree=zero 651 ekin=zero 652 do iatom=1,natom 653 do idir=1,3 654 ! Warning : the fixing of atomis is implemented in reduced 655 ! coordinates, so that this expression is wrong 656 if (iatfix(idir,iatom) == 0) then 657 ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2 658 ! Counts the degrees of freedom 659 nfree=nfree+one 660 end if 661 end do 662 end do 663 gnkt=nfree*ktemp 664 !Update the forces 665 glogs(1)=(two*ekin-gnkt)/qmass(1) 666 vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 667 do inos=1,nnos-1 668 alocal=exp(-dtion/eight*vlogs(nnos+1-inos)) 669 vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+& 670 & dtion/four*glogs(nnos-inos)*alocal 671 end do 672 !Update the particle velocities 673 alocal=exp(-dtion/two*vlogs(1)) 674 scale=scale*alocal 675 !Update the forces 676 glogs(1)=(scale*scale*two*ekin-gnkt)/qmass(1) 677 !Update the thermostat positions 678 do inos=1,nnos 679 xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two 680 end do 681 !Update the thermostat velocities 682 do inos=1,nnos-1 683 alocal=exp(-dtion/eight*vlogs(inos+1)) 684 vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal 685 glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1) 686 end do 687 vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four 688 vel(:,:)=vel(:,:)*scale 689 !Compute the ionic kinetic energy 690 ekin=zero 691 do iatom=1,natom 692 do idir=1,3 693 ! Warning : the fixing of atomis is implemented in reduced 694 ! coordinates, so that this expression is wrong 695 if (iatfix(idir,iatom) == 0) then 696 ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2 697 end if 698 end do 699 end do 700 !Compute the thermostat kinetic energy and add it to the ionic one 701 ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*nfree*ktemp 702 do inos=2,nnos 703 ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp 704 end do 705 mttk_vars%glogs(:)=glogs(:) 706 mttk_vars%vlogs(:)=vlogs(:) 707 mttk_vars%xlogs(:)=xlogs(:) 708 ABI_FREE(glogs) 709 ABI_FREE(vlogs) 710 ABI_FREE(xlogs) 711 !DEBUG 712 !write(std_out,*)'ekin added',half*qmass(1)*vlogs(1)**2,xlogs(1)*(nfree)*ktemp 713 !ENDEBUG 714 715 end subroutine isotemp
ABINIT/m_pred_isothermal [ Modules ]
NAME
m_pred_isothermal
FUNCTION
Ionmov predictors (13) Isothermal integrator 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 [[cite:Martyna1996]]
COPYRIGHT
Copyright (C) 1998-2024 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 .
SOURCE
20 #if defined HAVE_CONFIG_H 21 #include "config.h" 22 #endif 23 24 #include "abi_common.h" 25 26 module m_pred_isothermal 27 28 use defs_basis 29 use m_errors 30 use m_abicore 31 use m_abimover 32 use m_abihist 33 use m_linalg_interfaces 34 35 use m_numeric_tools, only : uniformrandom 36 use m_geometry, only : mkrdim, xcart2xred, xred2xcart, metric 37 38 implicit none 39 40 private
ABINIT/pred_isothermal [ Functions ]
NAME
pred_isothermal
FUNCTION
Ionmov predictors (13) Isothermal integrator IONMOV 13: Reversible integrator of Martyna at al. The equation of motion of the ions in contact with a thermostat and a barostat are solved with the algorithm proposed by Martyna, Tuckermann Tobias and Klein, Mol. Phys., 1996, p. 1117. [[cite:Martyna1996]] Related parameters : the time step (dtion), the initial temperature mdtemp(1), the final temperature mdtemp(2), the number of thermostats (nnos), and the masses of thermostats (qmass). If optcell=1 or 2, the mass of the barostat (bmass) must be given in addition. There are three sub cases according to the value of optcell optcell=0: isothermal optcell=1: homogeneous cell fluctuations optcell=2: full cell fluctuation in addition to temperature control.
INPUTS
ab_mover <type(abimover)> : Datatype with all the information needed by the preditor itime : Index of the present iteration ntime : Maximal number of iterations zDEBUG : if true print some debugging information
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
SOURCE
84 subroutine pred_isothermal(ab_mover,hist,itime,mttk_vars,ntime,zDEBUG,iexit) 85 86 implicit none 87 88 !Arguments ------------------------------------ 89 !scalars 90 type(abimover),intent(in) :: ab_mover 91 type(abihist),intent(inout) :: hist 92 type(mttk_type),intent(inout) :: mttk_vars 93 integer,intent(in) :: itime 94 integer,intent(in) :: ntime 95 integer,intent(in) :: iexit 96 logical,intent(in) :: zDEBUG 97 98 !Local variables------------------------------- 99 !scalars 100 integer :: ii,kk,iatom,idim,idum=5,ierr 101 integer,parameter :: lwork=8 102 real(dp) :: ucvol,ucvol0,ucvol_next,mttk_aloc,mttk_aloc2,mttk_bloc,ekin 103 real(dp) :: massvol=0 104 real(dp),parameter :: esh2=one/six,esh4=esh2/20._dp,esh6=esh4/42._dp 105 real(dp),parameter :: esh8=esh6/72._dp,nosetol=tol10,v2tol=tol8 106 real(dp) :: etotal,rescale_vel,polysh,s1,s2,sigma2,v2gauss,vtest 107 real(dp),save :: ktemp,vlogv 108 character(len=5000) :: message 109 !arrays 110 real(dp),allocatable,save :: fcart_m(:,:),vel_nexthalf(:,:) 111 112 real(dp) :: mttk_alc(3),mttk_alc2(3),mttk_blc(3),mttk_psh(3) 113 real(dp) :: mttk_tv(3,3),mttk_vt(3,3),mttk_ubox(3,3) 114 real(dp) :: mttk_uu(3),mttk_uv(3),mttk_veig(3) 115 real(dp) :: acell(3),acell0(3),acell_next(3) 116 real(dp) :: rprimd(3,3),rprimd0(3,3),rprim(3,3),rprimd_next(3,3),rprim_next(3,3) 117 real(dp) :: gprimd(3,3) 118 real(dp) :: gmet(3,3) 119 real(dp) :: rmet(3,3) 120 real(dp) :: fcart(3,ab_mover%natom) 121 real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom) 122 real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom) 123 real(dp) :: vel(3,ab_mover%natom) 124 real(dp) :: strten(6),work(lwork) 125 126 !*************************************************************************** 127 !Beginning of executable session 128 !*************************************************************************** 129 130 if(iexit/=0)then 131 if (allocated(fcart_m)) then 132 ABI_FREE(fcart_m) 133 end if 134 if (allocated(vel_nexthalf)) then 135 ABI_FREE(vel_nexthalf) 136 end if 137 return 138 end if 139 140 !write(std_out,*) 'isothermal 01' 141 !########################################################## 142 !### 01. Debugging and Verbose 143 144 if(zDEBUG)then 145 write(std_out,'(a,3a,41a,36a)') ch10,('-',kk=1,3),& 146 & 'Debugging and Verbose for pred_isothermal',('-',kk=1,36) 147 write(std_out,*) 'ionmov: ',13 148 write(std_out,*) 'itime: ',itime 149 end if 150 151 !write(std_out,*) 'isothermal 01' 152 !########################################################## 153 !### 01. Allocate the vectors vin, vout and hessian matrix 154 !### These arrays could be allocated from a previus 155 !### dataset that exit before itime==ntime 156 157 if(itime==1)then 158 if (allocated(fcart_m)) then 159 ABI_FREE(fcart_m) 160 end if 161 if (allocated(vel_nexthalf)) then 162 ABI_FREE(vel_nexthalf) 163 end if 164 end if 165 166 if (.not.allocated(fcart_m)) then 167 ABI_MALLOC(fcart_m,(3,ab_mover%natom)) 168 end if 169 if (.not.allocated(vel_nexthalf)) then 170 ABI_MALLOC(vel_nexthalf,(3,ab_mover%natom)) 171 end if 172 173 !write(std_out,*) 'isothermal 02' 174 !########################################################## 175 !### 02. Obtain the present values from the history 176 177 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 178 179 fcart(:,:)=hist%fcart(:,:,hist%ihist) 180 strten(:) =hist%strten(:,hist%ihist) 181 vel(:,:) =hist%vel(:,:,hist%ihist) 182 etotal =hist%etot(hist%ihist) 183 184 do ii=1,3 185 rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3) 186 end do 187 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 188 189 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 190 191 if(zDEBUG)then 192 write (std_out,*) 'fcart:' 193 do kk=1,ab_mover%natom 194 write (std_out,*) fcart(:,kk) 195 end do 196 write (std_out,*) 'vel:' 197 do kk=1,ab_mover%natom 198 write (std_out,*) vel(:,kk) 199 end do 200 write (std_out,*) 'strten:' 201 write (std_out,*) strten(1:3),ch10,strten(4:6) 202 write (std_out,*) 'etotal:' 203 write (std_out,*) etotal 204 end if 205 206 !Save initial values 207 acell0(:)=acell(:) 208 rprimd0(:,:)=rprimd(:,:) 209 ucvol0=ucvol 210 211 !write(std_out,*) 'isothermal 03' 212 !########################################################## 213 !### 05. Seconde half velocity step 214 215 if (itime>1) then 216 217 ! Next Half velocity step 218 do idim=1,3 219 fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:) 220 end do 221 vel(:,:)=vel_nexthalf(:,:)+ab_mover%dtion/two*fcart_m(:,:) 222 223 if (ab_mover%optcell==0) then 224 ! Update Thermostat variables and velocity 225 call isotemp(ab_mover%amass,ab_mover%dtion,ekin,ab_mover%iatfix,& 226 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,vel) 227 else if (ab_mover%optcell==1) then 228 ! Update Thermostat variables and velocity 229 call isopress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,& 230 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,& 231 & strten,ab_mover%strtarget,ucvol,vel,vlogv) 232 else if (ab_mover%optcell==2) then 233 ! Next half step for extended variables 234 call isostress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,& 235 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,& 236 & ab_mover%qmass,strten,ab_mover%strtarget,ucvol,vel) 237 end if 238 239 if(itime==2) massvol=ekin+etotal 240 241 if (ab_mover%optcell==2) then 242 ! Evolution of cell and volume 243 acell_next(:)=acell(:) 244 ucvol_next=ucvol 245 end if 246 247 call mkrdim(acell,rprim,rprimd) 248 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 249 250 if(zDEBUG)then 251 write(std_out,*) 'Second half velocity step' 252 write(std_out,*) 'Cell parameters:' 253 write(std_out,*) 'rprimd:' 254 do kk=1,3 255 write(std_out,*) rprimd(:,kk) 256 end do 257 write(std_out,*) 'rprim:' 258 do kk=1,3 259 write(std_out,*) rprim(:,kk) 260 end do 261 write(std_out,*) 'acell:' 262 write(std_out,*) acell(:) 263 write(std_out,*) 'Conserved energy:',(ekin+etotal)-massvol,ekin,etotal 264 write(std_out,*) 'Volume of unitqry cell (ucvol):',ucvol 265 end if 266 267 end if ! if (itime>1) 268 269 !write(std_out,*) 'isothermal 04' 270 !########################################################## 271 !### 03. Compute the next values 272 273 !The temperature is linear between initial and final values 274 !It is here converted from Kelvin to Hartree (kb_HaK) 275 ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK 276 277 if(zDEBUG)then 278 write(std_out,*) 'Temperature in Kelvin (ktemp):',ktemp 279 write(std_out,*) 'Initial temp (mdtemp(1)):',ab_mover%mdtemp(1) 280 write(std_out,*) 'Final temp (mdtemp(2)):',ab_mover%mdtemp(2) 281 write(std_out,*) 'Delay for atom permutation (delayperm)',ab_mover%delayperm 282 write(std_out,*) 'dtion',ab_mover%dtion 283 write(std_out,*) 'nnos:', ab_mover%nnos 284 write(std_out,*) 'qmass', ab_mover%qmass(:) 285 write(std_out,*) 'bmass',ab_mover%bmass 286 end if 287 288 if(itime==1) then 289 mttk_vars%glogs(:)=zero; mttk_vars%vlogs(:)=zero; mttk_vars%xlogs(:)=zero 290 mttk_vars%vboxg(:,:)=zero 291 vlogv=zero 292 ! v2gauss is twice the kinetic energy 293 v2gauss=0.0_dp 294 do iatom=1,ab_mover%natom 295 do idim=1,3 296 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 297 end do 298 end do 299 300 ! If there is no kinetic energy 301 if (v2gauss<=v2tol.and.itime==1) then 302 ! Maxwell-Boltzman distribution 303 v2gauss=zero 304 vtest=zero 305 do iatom=1,ab_mover%natom 306 do idim=1,3 307 vel(idim,iatom)=sqrt(kb_HaK*ab_mover%mdtemp(1)/ab_mover%amass(iatom))*cos(two_pi*uniformrandom(idum)) 308 vel(idim,iatom)=vel(idim,iatom)*sqrt(-2._dp*log(uniformrandom(idum))) 309 end do 310 end do 311 312 ! Get rid of center-of-mass velocity 313 s1=sum(ab_mover%amass(:)) 314 do idim=1,3 315 s2=sum(ab_mover%amass(:)*vel(idim,:)) 316 vel(idim,:)=vel(idim,:)-s2/s1 317 end do 318 319 ! Recompute v2gauss 320 do iatom=1,ab_mover%natom 321 do idim=1,3 322 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 323 vtest=vtest+vel(idim,iatom)/(3._dp*ab_mover%natom) 324 end do 325 end do 326 327 ! Now rescale the velocities to give the exact temperature 328 rescale_vel=sqrt(3._dp*ab_mover%natom*kb_HaK*ab_mover%mdtemp(1)/v2gauss) 329 vel(:,:)=vel(:,:)*rescale_vel 330 331 ! Recompute v2gauss with the rescaled velocities 332 v2gauss=zero 333 do iatom=1,ab_mover%natom 334 do idim=1,3 335 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 336 end do 337 end do 338 339 ! Compute the variance and print 340 sigma2=(v2gauss/(3._dp*ab_mover%natom)-ab_mover%amass(1)*vtest**2)/kb_HaK 341 342 if (zDEBUG)then 343 write(message, '(a)' )& 344 & ' Rescaling or initializing velocities to initial temperature' 345 call wrtout(std_out,message,'COLL') 346 write(message, '(a,d12.5,a,D12.5)' )& 347 & ' --- Scaling factor :',rescale_vel,' Asked T (K) ',ab_mover%mdtemp(1) 348 call wrtout(std_out,message,'COLL') 349 write(message, '(a,d12.5,a,D12.5)' )& 350 & ' --- Effective temperature',v2gauss/(3*ab_mover%natom*kb_HaK),' From variance', sigma2 351 call wrtout(std_out,message,'COLL') 352 end if 353 354 end if !(v2gauss<=v2tol.and.itime==1) 355 end if !(itime==1) 356 357 !XG070613 : Do not take away the following line , seems needed for the pathscale compiler 358 359 if (zDEBUG) write(std_out,*) 'vboxg',mttk_vars%vboxg(:,:) 360 361 362 !write(std_out,*) 'isothermal 05' 363 !########################################################## 364 !### 03. First half velocity step 365 366 !write(std_out,*) 'FIRST HALF VELOCITY STEP',ucvol 367 !write(std_out,*) 'OPTCELL option selected:',ab_mover%optcell 368 !write(std_out,*) 'RPRIMD' 369 !do kk=1,3 370 !write(std_out,*) rprimd(:,kk) 371 !end do 372 !write(std_out,*) 'RPRIM' 373 !do kk=1,3 374 !write(std_out,*) rprim(:,kk) 375 !end do 376 !write(std_out,*) 'ACELL' 377 !write(std_out,*) acell(:) 378 379 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 380 !%%% BEGIN sub case optcell=0 Isothermal Ensemble 381 !%%% 382 if(ab_mover%optcell==0) then 383 ! There is no evolution of cell 384 acell_next(:)=acell(:) 385 ucvol_next=ucvol 386 rprim_next(:,:)=rprim(:,:) 387 rprimd_next(:,:)=rprimd(:,:) 388 ! Update Thermostat variables and scale velocitie 389 call isotemp(ab_mover%amass,ab_mover%dtion,ekin,ab_mover%iatfix,& 390 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,vel) 391 392 ! Half velocity step 393 do idim=1,3 394 fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:) 395 end do 396 vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:) 397 ! New positions 398 ! Convert input xred (reduced coordinates) to xcart (cartesian) 399 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 400 xcart_next(:,:)=xcart(:,:)+vel_nexthalf(:,:)*ab_mover%dtion 401 ! Convert back to xred (reduced coordinates) 402 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 403 ! %%% 404 ! %%% END sub case optcell=0 Isothermal Ensemble 405 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 406 407 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 408 ! %%% BEGIN sub case optcell=1 Isothermal-Isenthalpic 409 ! %%% Ensemble (homogeneous cell deformation) 410 ! %%% 411 else if (ab_mover%optcell==1) then 412 ! Only homogeneous evolution of cell 413 ! Evolution of cell we keep rprim constant 414 rprim_next(:,:)=rprim(:,:) 415 ! Update Thermostat variables and velocity 416 call isopress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,& 417 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,& 418 & strten,ab_mover%strtarget,ucvol,vel,vlogv) 419 420 ! Half velocity step 421 do idim=1,3 422 fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:) 423 end do 424 vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:) 425 ! New positions 426 mttk_aloc=exp(ab_mover%dtion/two*vlogv) 427 mttk_aloc2=(vlogv*ab_mover%dtion/two)**2 428 polysh=(((esh8*mttk_aloc2+esh6)*mttk_aloc2+esh4)*mttk_aloc2+esh2)*mttk_aloc2+one 429 mttk_bloc=mttk_aloc*polysh*ab_mover%dtion 430 ! Convert input xred (reduced coordinates) to xcart (cartesian) 431 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 432 xcart_next(:,:)=xcart(:,:)*mttk_aloc**2+vel_nexthalf(:,:)*mttk_bloc 433 ! Update the volume and related quantities 434 acell_next(:)=acell(:)*exp(ab_mover%dtion*vlogv) 435 ! ucvol=ucvol*exp(ab_mover%dtion*vlogv) 436 call mkrdim(acell_next,rprim,rprimd_next) 437 call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol_next) 438 ! Convert back to xred (reduced coordinates) 439 call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next) 440 ! Computation of the forces for the new positions 441 ! Compute DFT forces (big loop) 442 443 ! COMMENTED 444 ! This should be in mover.F90 445 446 ! ! If metric has changed since the initialization, update the Ylm's 447 ! if (ab_mover%optcell/=0.and.psps%useylm==1.and.itime>1)then 448 ! option=0;if (ab_mover%iscf>0) option=1 449 ! call initylmg(gprimd,kg,ab_mover%kptns,ab_mover%mkmem,mpi_enreg,psps%mpsang,ab_mover%mpw,ab_mover%nband,ab_mover%nkpt,& 450 ! & npwarr,ab_mover%nsppol,option,rprimd_next,ylm,ylmgr) 451 ! end if 452 453 454 ! %%% 455 ! %%% END sub case optcell=1 Isothermal-Isenthalpic 456 ! %%% Ensemble (homogeneous cell deformation) 457 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 458 459 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 460 ! %%% BEGIN sub case optcell=2 Isothermal-Isenthalpic 461 ! %%% Ensemble (full cell deformation) 462 ! %%% 463 else if (ab_mover%optcell==2) then 464 acell_next=acell 465 ! Fisrt half step for extended variables 466 call isostress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,& 467 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,& 468 & ab_mover%qmass,strten,ab_mover%strtarget,ucvol,vel) 469 ! Half velocity step 470 do idim=1,3 471 fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:) 472 end do 473 vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:) 474 ! Convert input xred (reduced coordinates) to xcart (cartesian) 475 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 476 ! New positions 477 mttk_vt(:,:)=mttk_vars%vboxg(:,:) 478 call dsyev('V','U',3,mttk_vt,3,mttk_veig,work,lwork,ierr) 479 mttk_tv(:,:)=transpose(mttk_vt) 480 mttk_alc(:)=exp(ab_mover%dtion/two*mttk_veig(:)) 481 mttk_alc2(:)=(mttk_veig(:)*ab_mover%dtion/two)**2 482 mttk_psh(:)=(((esh8*mttk_alc2(:)+esh6)*mttk_alc2(:)+esh4)*mttk_alc2(:)+esh2)*mttk_alc2(:)+one 483 mttk_blc(:)=mttk_alc(:)*mttk_psh(:)*ab_mover%dtion 484 ! Update the positions 485 do iatom=1,ab_mover%natom 486 mttk_uu(:)=matmul(mttk_tv,xcart(:,iatom)) 487 mttk_uv(:)=matmul(mttk_tv,vel_nexthalf(:,iatom)) 488 mttk_uu(:)=mttk_uu(:)*mttk_alc(:)**2+mttk_uv(:)*mttk_blc(:) 489 xcart_next(:,iatom)=matmul(mttk_vt,mttk_uu) 490 end do 491 ! Update the box (rprimd and rprim) 492 mttk_ubox(:,:)=matmul(mttk_tv,rprimd) 493 do idim=1,3 494 mttk_ubox(:,idim)=mttk_ubox(:,idim)*mttk_alc(:)**2 495 end do 496 rprimd_next(:,:)=matmul(mttk_vt,mttk_ubox) 497 do idim=1,3 498 rprim_next(idim,:)=rprimd_next(idim,:)/acell(:) 499 end do 500 ! Update the volume 501 call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol) 502 ! Convert back to xred (reduced coordinates) 503 call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next) 504 ! Computation of the forces for the new positions 505 506 ! COMMENTED 507 ! This should be in mover.F90 508 509 ! ! If metric has changed since the initialization, update the Ylm's 510 ! if (ab_mover%optcell/=0.and.psps%useylm==1.and.itime>1)then 511 ! option=0;if (ab_mover%iscf>0) option=1 512 ! call initylmg(gprimd,kg,ab_mover%kptns,ab_mover%mkmem,mpi_enreg,psps%mpsang,ab_mover%mpw,ab_mover%nband,ab_mover%nkpt,& 513 ! & npwarr,ab_mover%nsppol,option,rprimd_next,ylm,ylmgr) 514 ! end if 515 516 ! %%% 517 ! %%% END sub case optcell=2 Isothermal-Isenthalpic 518 ! %%% Ensemble (full cell deformation) 519 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 520 else 521 write(message, '(a,i12,a,a)' )& 522 & ' Disallowed value for optcell=',ab_mover%optcell,ch10,& 523 & ' Allowed values with ionmov==13 : 0 to 2.' 524 ABI_BUG(message) 525 end if 526 527 !write(std_out,*) 'OLD PARAMETERS' 528 !write(std_out,*) 'RPRIMD' 529 !do kk=1,3 530 !write(std_out,*) rprimd(:,kk) 531 !end do 532 !write(std_out,*) 'RPRIM' 533 !do kk=1,3 534 !write(std_out,*) rprim(:,kk) 535 !end do 536 !write(std_out,*) 'ACELL' 537 !write(std_out,*) acell(:) 538 539 !write(std_out,*) 'NEXT PARAMETERS' 540 !write(std_out,*) 'RPRIMD' 541 !do kk=1,3 542 !write(std_out,*) rprimd_next(:,kk) 543 !end do 544 !write(std_out,*) 'RPRIM' 545 !do kk=1,3 546 !write(std_out,*) rprim_next(:,kk) 547 !end do 548 !write(std_out,*) 'ACELL' 549 !write(std_out,*) acell_next(:) 550 551 552 !Those are the values store into the history 553 rprim=rprim_next 554 rprimd=rprimd_next 555 xred=xred_next 556 xcart=xcart_next 557 acell=acell_next 558 559 !write(std_out,*) 'isothermal 06' 560 !########################################################## 561 !### 06. Update the history with the prediction 562 563 !Increase indexes 564 hist%ihist = abihist_findIndex(hist,+1) 565 566 !Fill the history with the variables 567 !xred, acell, rprimd, vel 568 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 569 hist%vel(:,:,hist%ihist)=vel(:,:) 570 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 571 572 if(zDEBUG)then 573 write (std_out,*) 'fcart:' 574 do kk=1,ab_mover%natom 575 write (std_out,*) fcart(:,kk) 576 end do 577 write (std_out,*) 'vel:' 578 do kk=1,ab_mover%natom 579 write (std_out,*) vel(:,kk) 580 end do 581 write (std_out,*) 'strten:' 582 write (std_out,*) strten(1:3),ch10,strten(4:6) 583 write (std_out,*) 'etotal:' 584 write (std_out,*) etotal 585 end if 586 587 end subroutine pred_isothermal