TABLE OF CONTENTS


ABINIT/isopress [ Functions ]

[ Top ] [ 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

PARENTS

      pred_isothermal

CHILDREN

      dsyev

SOURCE

799  subroutine isopress(amass,bmass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,qmass,&
800    & strten,strtarget,ucvol,vel,vlogv)
801 
802 
803 !This section has been created automatically by the script Abilint (TD).
804 !Do not modify the following lines by hand.
805 #undef ABI_FUNC
806 #define ABI_FUNC 'isopress'
807 !End of the abilint section
808 
809  implicit none
810 
811 !Arguments ------------------------------------
812 !scalars
813  integer,intent(in) :: nnos,natom
814  real(dp),intent(in) :: dtion,ktemp,ucvol,bmass
815  real(dp),intent(inout) :: vlogv
816  real(dp),intent(out) :: ekin
817  type(mttk_type) :: mttk_vars
818 !arrays
819  real(dp),intent(in) :: amass(natom),strtarget(6),strten(6)
820  real(dp),intent(inout) :: vel(3,natom)
821  real(dp),intent(in) :: qmass(:)
822  integer,intent(in) :: iatfix(:,:)
823 
824 !Local variables ------------------------------
825 !scalars
826  integer :: iatom,idir,inos
827  real(dp) :: alocal,glogv,gn1kt,gnkt,nfree,odnf,press,prtarget,scale
828  logical :: DEBUG=.FALSE.
829 !character(len=500) :: message
830 !arrays
831  real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:)
832 
833 !***************************************************************************
834 !Beginning of executable session
835 !***************************************************************************
836 
837  if(DEBUG) then
838    write(std_out,*) ch10,'***',ch10
839    write(std_out,*)' isopress : enter '
840    write(std_out,*)' strtarget=',strtarget
841    write(std_out,*)' bmass=',bmass
842    write(std_out,*)' dtion=',dtion
843    write(std_out,*)' ktemp=',ktemp
844    write(std_out,*)' natom=',natom
845    write(std_out,*)' nnos=',nnos
846    write(std_out,*)' strten=',strten
847    write(std_out,*)' ucvol=',ucvol
848  end if
849 
850  ABI_ALLOCATE(glogs,(nnos))
851  ABI_ALLOCATE(vlogs,(nnos))
852  ABI_ALLOCATE(xlogs,(nnos))
853  glogs(:)=mttk_vars%glogs(:)
854  vlogs(:)=mttk_vars%vlogs(:)
855  xlogs(:)=mttk_vars%xlogs(:)
856  glogv   =mttk_vars%glogv
857  scale=one
858 !Compute the ionic kinetic energy
859  nfree=zero
860  ekin=zero
861  do iatom=1,natom
862    do idir=1,3
863 !    Warning : the fixing of atomis is implemented in reduced
864 !    coordinates, so that this expression is wrong
865      if (iatfix(idir,iatom) == 0) then
866        ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2
867 !      Counts the degrees of freedom
868        nfree=nfree+one
869      end if
870    end do
871  end do
872  prtarget=-(strtarget(1)+strtarget(2)+strtarget(3))/three
873  press=-(strten(1)+strten(2)+strten(3))/three
874  gnkt=nfree*ktemp
875  gn1kt=(nfree+one)*ktemp
876  odnf=one+three/nfree
877 !Update the forces
878  glogs(1)=(two*ekin+bmass*vlogv*vlogv-gn1kt)/qmass(1)
879  glogv=(odnf*two*ekin+three*(press-prtarget)*ucvol)/bmass
880 !Update thermostat velocity
881  vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
882  do inos=1,nnos-1
883    alocal=exp(-dtion/eight*vlogs(nnos+1-inos))
884    vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+&
885 &   dtion/four*glogs(nnos-inos)*alocal
886  end do
887 !Update dLog(V)/dt
888  alocal=exp(-dtion/eight*vlogs(1))
889  vlogv=vlogv*alocal**2+dtion/four*glogv*alocal
890 !Update the particle velocities
891  alocal=exp(-dtion/two*(vlogs(1)+odnf*vlogv))
892  scale=scale*alocal
893  ekin=ekin*alocal**2
894  glogv=(odnf*two*ekin+three*(press-prtarget)*ucvol)/bmass
895 !Update the thermostat positions
896  do inos=1,nnos
897    xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two
898  end do
899 !Update dLog(V)/dt
900  alocal=exp(-dtion/eight*vlogs(1))
901  vlogv=vlogv*alocal**2+dtion/four*glogv*alocal
902 !Update the forces
903  glogs(1)=(two*ekin+bmass*vlogv*vlogv-gn1kt)/qmass(1)
904 !Update the thermostat velocities
905  do inos=1,nnos-1
906    alocal=exp(-dtion/eight*vlogs(inos+1))
907    vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal
908    glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1)
909  end do
910  vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
911  vel(:,:)=vel(:,:)*scale
912 !Compute the ionic kinetic energy
913  ekin=zero
914  do iatom=1,natom
915    do idir=1,3
916 !    Warning : the fixing of atomis is implemented in reduced
917 !    coordinates, so that this expression is wrong
918      if (iatfix(idir,iatom) == 0) then
919        ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2
920      end if
921    end do
922  end do
923 !Compute the thermostat kinetic energy and add it to the ionic one
924 !First thermostat
925  ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*(nfree+one)*ktemp
926 !Other thermostats
927  do inos=2,nnos
928    ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp
929  end do
930 !Barostat
931  ekin=ekin+half*bmass*vlogv**2+prtarget*ucvol
932  ABI_DEALLOCATE(glogs)
933  ABI_DEALLOCATE(vlogs)
934  ABI_DEALLOCATE(xlogs)
935 
936 !DEBUG
937 !write(std_out,*) 'EKIN',ekin
938 !write(std_out,*) 'VLOGV',vlogv
939 !write(std_out,*)'ekin added T',half*qmass(:)*vlogs(:)**2,xlogs(:)*(nfree)*ktemp
940 !write(std_out,*)'ekin added P',half*bmass*vlogv**2,prtarget*ucvol
941 !write(std_out,*)'ekin last',ekin
942 !ENDDEBUG
943 end subroutine isopress

ABINIT/isostress [ Functions ]

[ Top ] [ 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

PARENTS

      pred_isothermal

CHILDREN

      dsyev

SOURCE

 980  subroutine isostress(amass,bmass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,&
 981    & qmass,strten,strtarget,ucvol,vel)
 982 
 983  use m_linalg_interfaces
 984 
 985 !This section has been created automatically by the script Abilint (TD).
 986 !Do not modify the following lines by hand.
 987 #undef ABI_FUNC
 988 #define ABI_FUNC 'isostress'
 989 !End of the abilint section
 990 
 991  implicit none
 992 
 993 !Arguments ------------------------------------
 994 !scalars
 995  integer,intent(in) :: natom,nnos
 996  real(dp),intent(in) :: dtion,ktemp,ucvol,bmass
 997  real(dp),intent(out) :: ekin
 998  type(mttk_type) :: mttk_vars
 999 !arrays
1000  real(dp),intent(in) :: amass(natom),strtarget(6),strten(6),qmass(:)
1001  real(dp),intent(inout) :: vel(3,natom)
1002  integer, intent(in) :: iatfix(:,:)
1003 
1004 !Local variables ------------------------------
1005 !scalars
1006  integer,parameter :: lwork=8
1007  integer :: iatom,idir,info,inos,jdir
1008  real(dp) :: akinb,alocal,gn1kt,gnd2kt,nfree,odnf,trvg !,gnkt,scale
1009  logical  :: DEBUG=.FALSE.
1010 !character(len=500) :: message
1011 !arrays
1012  real(dp) :: akin(3,3),expdiag(3),gboxg(3,3),identity(3,3),press(3,3)
1013  real(dp) :: prtarget(3,3),tvtemp(3,3),uv(3),vboxg(3,3),veig(3),vtemp(3,3)
1014  real(dp) :: work(lwork)
1015  real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:)
1016 
1017 !***************************************************************************
1018 !Beginning of executable session
1019 !***************************************************************************
1020 
1021  if(DEBUG) then
1022    write(std_out,*) ch10,'***',ch10
1023    write(std_out,*)' isostress : enter '
1024    write(std_out,*)' strtarget=',strtarget
1025    write(std_out,*)' bmass=',bmass
1026    write(std_out,*)' dtion=',dtion
1027    write(std_out,*)' ktemp=',ktemp
1028    write(std_out,*)' natom=',natom
1029    write(std_out,*)' nnos=',nnos
1030    write(std_out,*)' strten=',strten
1031    write(std_out,*)' ucvol=',ucvol
1032  end if
1033 
1034  ABI_ALLOCATE(glogs,(nnos))
1035  ABI_ALLOCATE(vlogs,(nnos))
1036  ABI_ALLOCATE(xlogs,(nnos))
1037  glogs(:)=mttk_vars%glogs(:)
1038  vlogs(:)=mttk_vars%vlogs(:)
1039  xlogs(:)=mttk_vars%xlogs(:)
1040  vboxg(:,:)=mttk_vars%vboxg(:,:)
1041  identity(:,:)=zero
1042  do idir=1,3
1043    identity(idir,idir)=one
1044  end do
1045 
1046 !write(std_out,*) 'isostress 02'
1047 !##########################################################
1048 !### 03. Compute the ionic kinetic energy
1049 
1050  nfree=zero
1051  ekin=zero
1052  do iatom=1,natom
1053    do idir=1,3
1054 !    Warning : the fixing of atomis is implemented in reduced
1055 !    coordinates, so that this exprtargetion is wrong
1056      if (iatfix(idir,iatom) == 0) then
1057        ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2
1058 !      Counts the degrees of freedom
1059        nfree=nfree+one
1060      end if
1061    end do
1062  end do
1063 
1064  gn1kt=(nfree+one)*ktemp
1065  gnd2kt=(nfree+9)*ktemp
1066  odnf=one+three/nfree
1067  akin(:,:)=zero
1068  do iatom=1,natom
1069    do idir=1,3
1070      do jdir=1,3
1071 !      Warning : the fixing of atomis is implemented in reduced
1072 !      coordinates, so that this expression is wrong
1073        akin(idir,jdir)=akin(idir,jdir)+0.5d0*amass(iatom)*vel(idir,iatom)*vel(jdir,iatom)
1074      end do
1075    end do
1076  end do
1077  akinb=zero
1078  do idir=1,3
1079    do jdir=1,3
1080      akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2
1081    end do
1082  end do
1083 !Compute the pressure: from Voigt to tensor notation+kinetic energy
1084  do idir=1,3
1085    press(idir,idir)=-strten(idir)
1086    prtarget(idir,idir)=-strtarget(idir)
1087  end do
1088  press(3,2)=-strten(4); press(1,3)=-strten(5); press(2,1)=-strten(6)
1089  prtarget(3,2)=-strtarget(4); prtarget(1,3)=-strtarget(5); prtarget(2,1)=-strtarget(6)
1090  press(2,3)=press(3,2); press(3,1)=press(1,3); press(1,2)=press(2,1)
1091  prtarget(2,3)=prtarget(3,2); prtarget(3,1)=prtarget(1,3); prtarget(1,2)=prtarget(2,1)
1092 !Update the forces
1093  glogs(1)=(two*ekin+two*akinb-gnd2kt)/qmass(1)
1094  gboxg(:,:)=(two*ekin/nfree*identity(:,:)+two*akin(:,:)+(press(:,:)-prtarget(:,:))*ucvol)/bmass
1095 !Update thermostat velocity
1096  if (nnos > 0) vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
1097  do inos=1,nnos-1
1098    alocal=exp(-dtion/eight*vlogs(nnos+1-inos))
1099    vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+&
1100 &   dtion/four*glogs(nnos-inos)*alocal
1101  end do
1102 !Update box velocity
1103  alocal=exp(-dtion/eight*vlogs(1))
1104 
1105  if(DEBUG) then
1106    write(std_out,*)' gboxg(:,:)=',gboxg(:,:)
1107    write(std_out,*)' vboxg(:,:)=',vboxg(:,:)
1108    write(std_out,*)' alocal=',alocal
1109  end if
1110 
1111  vboxg(:,:)=vboxg(:,:)*alocal**2+dtion/four*gboxg(:,:)*alocal
1112 !Update the thermostat positions
1113  do inos=1,nnos
1114    xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two
1115  end do
1116 !Update the particle velocities
1117  trvg=(vboxg(1,1)+vboxg(2,2)+vboxg(3,3))/nfree
1118  vtemp(:,:)=vboxg(:,:)+(trvg+vlogs(1))*identity(:,:)
1119  call dsyev('V','U',3,vtemp,3,veig,work,lwork,info)
1120 !On exit, we have vtemp=U such that tU vtemp U = veig
1121  tvtemp(:,:)=transpose(vtemp)
1122 
1123  if(DEBUG) then
1124    write(std_out,*)' vboxg(:,:)=',vboxg(:,:)
1125    write(std_out,*)' vtemp(:,:)=',vtemp(:,:)
1126    write(std_out,*)' veig(:)=',veig(:)
1127  end if
1128 
1129  expdiag(1)=exp(-veig(1)*dtion/two)
1130  expdiag(2)=exp(-veig(2)*dtion/two)
1131  expdiag(3)=exp(-veig(3)*dtion/two)
1132  if(DEBUG) then
1133    write(std_out,*)' isostress : expdiag(:)=',expdiag(:)  ! Do not remove this line : seems to be needed for g95 compilo
1134  end if
1135  do iatom=1,natom
1136    uv(:)=matmul(tvtemp,vel(:,iatom))
1137    uv(:)=uv(:)*expdiag(:)
1138    vel(:,iatom)=matmul(vtemp,uv)
1139  end do
1140 !Compute the ionic kinetic energy
1141  nfree=zero
1142  ekin=zero
1143  do iatom=1,natom
1144    do idir=1,3
1145 !    Warning : the fixing of atomis is implemented in reduced
1146 !    coordinates, so that this expression is wrong
1147      if (iatfix(idir,iatom) == 0) then
1148        ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2
1149        if(DEBUG) then
1150          write(std_out,*)'kin',iatom,ekin,vel(idir,iatom)
1151        end if
1152 !      Counts the degrees of freedom
1153        nfree=nfree+one
1154      end if
1155    end do
1156  end do
1157  gn1kt=(nfree+one)*ktemp
1158  gnd2kt=(nfree+9)*ktemp
1159  odnf=one+three/nfree
1160  akin(:,:)=zero
1161  do iatom=1,natom
1162    do idir=1,3
1163      do jdir=1,3
1164 !      Warning : the fixing of atomis is implemented in reduced
1165 !      coordinates, so that this expression is wrong
1166        akin(idir,jdir)=akin(idir,jdir)+0.5d0*amass(iatom)*vel(idir,iatom)*vel(jdir,iatom)
1167      end do
1168    end do
1169  end do
1170  gboxg(:,:)=(two*ekin/nfree*identity(:,:)+two*akin(:,:)+(press(:,:)-prtarget(:,:))*ucvol)/bmass
1171 !Update box velocity
1172  alocal=exp(-dtion/eight*vlogs(1))
1173  vboxg(:,:)=vboxg(:,:)*alocal**2+dtion/four*gboxg(:,:)*alocal
1174 !Compute the box kinetic energy
1175  akinb=zero
1176  do idir=1,3
1177    do jdir=1,3
1178      akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2
1179    end do
1180  end do
1181  glogs(1)=(two*ekin+two*akinb-gnd2kt)/qmass(1)
1182 !Update the thermostat velocities
1183  do inos=1,nnos-1
1184    alocal=exp(-dtion/eight*vlogs(inos+1))
1185    vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal
1186    glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1)
1187  end do
1188  if (nnos > 0) vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
1189 !Compute the ionic kinetic energy
1190  ekin=zero
1191  do iatom=1,natom
1192    do idir=1,3
1193 !    Warning : the fixing of atomis is implemented in reduced
1194 !    coordinates, so that this expression is wrong
1195      if (iatfix(idir,iatom) == 0) then
1196        ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2
1197      end if
1198    end do
1199  end do
1200 !Compute the thermostat kinetic energy and add it to the ionic one
1201 !First thermostat
1202  ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*(nfree+nine)*ktemp
1203 !Other thermostats
1204  do inos=2,nnos
1205    ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp
1206  end do
1207 !Barostat kinetic energy
1208  akinb=zero
1209  do idir=1,3
1210    do jdir=1,3
1211      akinb=akinb+0.5d0*bmass*vboxg(idir,jdir)**2
1212    end do
1213  end do
1214 !ekin is the invariant minus the potential energy
1215  ekin=ekin+akinb+prtarget(1,1)*ucvol
1216 
1217  mttk_vars%vboxg(:,:)=vboxg(:,:)
1218  mttk_vars%glogs(:)=glogs(:)
1219  mttk_vars%vlogs(:)=vlogs(:)
1220  mttk_vars%xlogs(:)=xlogs(:)
1221  ABI_DEALLOCATE(glogs)
1222  ABI_DEALLOCATE(vlogs)
1223  ABI_DEALLOCATE(xlogs)
1224 
1225  if(DEBUG) then
1226 !  write(std_out,*)'ekin added T',half*qmass(:)*vlogs(:)**2,xlogs(:)*(nfree)*ktemp
1227 !  write(std_out,*)'ekin added P',akinb,prtarget*ucvol
1228 !  write(std_out,*)'ekin last',ekin
1229    write(std_out,*) ch10,' exiting from isostress',ch10
1230  end if
1231 
1232 end subroutine isostress

ABINIT/isotemp [ Functions ]

[ Top ] [ 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

PARENTS

      pred_isothermal

CHILDREN

      dsyev

SOURCE

653 subroutine isotemp(amass,dtion,ekin,iatfix,ktemp,mttk_vars,natom,nnos,qmass,vel)
654 
655 
656 !This section has been created automatically by the script Abilint (TD).
657 !Do not modify the following lines by hand.
658 #undef ABI_FUNC
659 #define ABI_FUNC 'isotemp'
660 !End of the abilint section
661 
662  implicit none
663 
664 !Arguments ------------------------------------
665 !scalars
666  integer,intent(in) :: natom
667  integer,intent(in) :: nnos
668  real(dp),intent(in) :: dtion,ktemp
669  real(dp),intent(out) :: ekin
670  type(mttk_type) :: mttk_vars
671 !arrays
672  real(dp),intent(in) :: amass(natom)
673  real(dp),intent(inout) :: vel(3,natom)
674  real(dp),intent(in) :: qmass(:)
675  integer,intent(in) :: iatfix(:,:)
676 
677 !Local variables ------------------------------
678 !scalars
679  integer :: iatom,idir,inos
680  real(dp) :: alocal,gnkt,nfree,scale
681  !character(len=500) :: message
682 !arrays
683  real(dp),allocatable :: glogs(:),vlogs(:),xlogs(:)
684 
685 !***************************************************************************
686 !Beginning of executable session
687 !***************************************************************************
688 
689  ABI_ALLOCATE(glogs,(nnos))
690  ABI_ALLOCATE(vlogs,(nnos))
691  ABI_ALLOCATE(xlogs,(nnos))
692  glogs(:)=mttk_vars%glogs(:)
693  vlogs(:)=mttk_vars%vlogs(:)
694  xlogs(:)=mttk_vars%xlogs(:)
695  scale=one
696 !Compute the ionic kinetic energy
697  nfree=zero
698  ekin=zero
699  do iatom=1,natom
700    do idir=1,3
701 !    Warning : the fixing of atomis is implemented in reduced
702 !    coordinates, so that this expression is wrong
703      if (iatfix(idir,iatom) == 0) then
704        ekin=ekin+0.5d0*amass(iatom)*vel(idir,iatom)**2
705 !      Counts the degrees of freedom
706        nfree=nfree+one
707      end if
708    end do
709  end do
710  gnkt=nfree*ktemp
711 !Update the forces
712  glogs(1)=(two*ekin-gnkt)/qmass(1)
713  vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
714  do inos=1,nnos-1
715    alocal=exp(-dtion/eight*vlogs(nnos+1-inos))
716    vlogs(nnos-inos)=vlogs(nnos-inos)*alocal*alocal+&
717 &   dtion/four*glogs(nnos-inos)*alocal
718  end do
719 !Update the particle velocities
720  alocal=exp(-dtion/two*vlogs(1))
721  scale=scale*alocal
722 !Update the forces
723  glogs(1)=(scale*scale*two*ekin-gnkt)/qmass(1)
724 !Update the thermostat positions
725  do inos=1,nnos
726    xlogs(inos)=xlogs(inos)+vlogs(inos)*dtion/two
727  end do
728 !Update the thermostat velocities
729  do inos=1,nnos-1
730    alocal=exp(-dtion/eight*vlogs(inos+1))
731    vlogs(inos)=vlogs(inos)*alocal*alocal+dtion/four*glogs(inos)*alocal
732    glogs(inos+1)=(qmass(inos)*vlogs(inos)*vlogs(inos)-ktemp)/qmass(inos+1)
733  end do
734  vlogs(nnos)=vlogs(nnos)+glogs(nnos)*dtion/four
735  vel(:,:)=vel(:,:)*scale
736 !Compute the ionic kinetic energy
737  ekin=zero
738  do iatom=1,natom
739    do idir=1,3
740 !    Warning : the fixing of atomis is implemented in reduced
741 !    coordinates, so that this expression is wrong
742      if (iatfix(idir,iatom) == 0) then
743        ekin=ekin+half*amass(iatom)*vel(idir,iatom)**2
744      end if
745    end do
746  end do
747 !Compute the thermostat kinetic energy and add it to the ionic one
748  ekin=ekin+half*qmass(1)*vlogs(1)**2+xlogs(1)*nfree*ktemp
749  do inos=2,nnos
750    ekin=ekin+half*qmass(inos)*vlogs(inos)**2+xlogs(inos)*ktemp
751  end do
752  mttk_vars%glogs(:)=glogs(:)
753  mttk_vars%vlogs(:)=vlogs(:)
754  mttk_vars%xlogs(:)=xlogs(:)
755  ABI_DEALLOCATE(glogs)
756  ABI_DEALLOCATE(vlogs)
757  ABI_DEALLOCATE(xlogs)
758 !DEBUG
759 !write(std_out,*)'ekin added',half*qmass(1)*vlogs(1)**2,xlogs(1)*(nfree)*ktemp
760 !ENDEBUG
761 
762 end subroutine isotemp

ABINIT/m_pred_isothermal [ Modules ]

[ Top ] [ 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-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 .

PARENTS

CHILDREN

SOURCE

25 #if defined HAVE_CONFIG_H
26 #include "config.h"
27 #endif
28 
29 #include "abi_common.h"
30 
31 module m_pred_isothermal
32 
33  use defs_basis
34  use m_errors
35  use m_abicore
36  use m_abimover
37  use m_abihist
38  use m_linalg_interfaces
39 
40  use m_numeric_tools,  only : uniformrandom
41  use m_geometry,       only : mkrdim, xcart2xred, xred2xcart, metric
42 
43  implicit none
44 
45  private

ABINIT/pred_isothermal [ Functions ]

[ Top ] [ 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

PARENTS

      mover

CHILDREN

      dsyev,hist2var,isopress,isostress,isotemp,metric,mkrdim,var2hist,wrtout
      xcart2xred,xred2xcart

SOURCE

 96 subroutine pred_isothermal(ab_mover,hist,itime,mttk_vars,ntime,zDEBUG,iexit)
 97 
 98 
 99 !This section has been created automatically by the script Abilint (TD).
100 !Do not modify the following lines by hand.
101 #undef ABI_FUNC
102 #define ABI_FUNC 'pred_isothermal'
103 !End of the abilint section
104 
105  implicit none
106 
107 !Arguments ------------------------------------
108 !scalars
109  type(abimover),intent(in)       :: ab_mover
110  type(abihist),intent(inout) :: hist
111  type(mttk_type),intent(inout) :: mttk_vars
112  integer,intent(in) :: itime
113  integer,intent(in) :: ntime
114  integer,intent(in) :: iexit
115  logical,intent(in) :: zDEBUG
116 
117 !Local variables-------------------------------
118 !scalars
119  integer  :: ii,kk,iatom,idim,idum=5,ierr
120  integer,parameter :: lwork=8
121  real(dp) :: ucvol,ucvol0,ucvol_next,mttk_aloc,mttk_aloc2,mttk_bloc,ekin
122  real(dp) :: massvol=0
123  real(dp),parameter :: esh2=one/six,esh4=esh2/20._dp,esh6=esh4/42._dp
124  real(dp),parameter :: esh8=esh6/72._dp,nosetol=tol10,v2tol=tol8
125  real(dp) :: etotal,rescale_vel,polysh,s1,s2,sigma2,v2gauss,vtest
126  real(dp),save :: ktemp,vlogv
127  character(len=5000) :: message
128 !arrays
129  real(dp),allocatable,save :: fcart_m(:,:),vel_nexthalf(:,:)
130 
131  real(dp) :: mttk_alc(3),mttk_alc2(3),mttk_blc(3),mttk_psh(3)
132  real(dp) :: mttk_tv(3,3),mttk_vt(3,3),mttk_ubox(3,3)
133  real(dp) :: mttk_uu(3),mttk_uv(3),mttk_veig(3)
134  real(dp) :: acell(3),acell0(3),acell_next(3)
135  real(dp) :: rprimd(3,3),rprimd0(3,3),rprim(3,3),rprimd_next(3,3),rprim_next(3,3)
136  real(dp) :: gprimd(3,3)
137  real(dp) :: gmet(3,3)
138  real(dp) :: rmet(3,3)
139  real(dp) :: fcart(3,ab_mover%natom)
140 !real(dp) :: fred_corrected(3,ab_mover%natom)
141  real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
142  real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
143  real(dp) :: vel(3,ab_mover%natom)
144  real(dp) :: strten(6),work(lwork)
145 
146 !***************************************************************************
147 !Beginning of executable session
148 !***************************************************************************
149 
150  if(iexit/=0)then
151    if (allocated(fcart_m))       then
152      ABI_DEALLOCATE(fcart_m)
153    end if
154    if (allocated(vel_nexthalf))  then
155      ABI_DEALLOCATE(vel_nexthalf)
156    end if
157    return
158  end if
159 
160 !write(std_out,*) 'isothermal 01'
161 !##########################################################
162 !### 01. Debugging and Verbose
163 
164  if(zDEBUG)then
165    write(std_out,'(a,3a,41a,36a)') ch10,('-',kk=1,3),&
166 &   'Debugging and Verbose for pred_isothermal',('-',kk=1,36)
167    write(std_out,*) 'ionmov: ',13
168    write(std_out,*) 'itime:  ',itime
169  end if
170 
171 !write(std_out,*) 'isothermal 01'
172 !##########################################################
173 !### 01. Allocate the vectors vin, vout and hessian matrix
174 !###     These arrays could be allocated from a previus
175 !###     dataset that exit before itime==ntime
176 
177  if(itime==1)then
178    if (allocated(fcart_m))       then
179      ABI_DEALLOCATE(fcart_m)
180    end if
181    if (allocated(vel_nexthalf))  then
182      ABI_DEALLOCATE(vel_nexthalf)
183    end if
184  end if
185 
186  if (.not.allocated(fcart_m))       then
187    ABI_ALLOCATE(fcart_m,(3,ab_mover%natom))
188  end if
189  if (.not.allocated(vel_nexthalf))  then
190    ABI_ALLOCATE(vel_nexthalf,(3,ab_mover%natom))
191  end if
192 
193 !write(std_out,*) 'isothermal 02'
194 !##########################################################
195 !### 02. Obtain the present values from the history
196 
197  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
198 
199  fcart(:,:)=hist%fcart(:,:,hist%ihist)
200  strten(:) =hist%strten(:,hist%ihist)
201  vel(:,:)  =hist%vel(:,:,hist%ihist)
202  etotal    =hist%etot(hist%ihist)
203 
204  do ii=1,3
205    rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
206  end do
207  call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
208 
209  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
210 
211  if(zDEBUG)then
212    write (std_out,*) 'fcart:'
213    do kk=1,ab_mover%natom
214      write (std_out,*) fcart(:,kk)
215    end do
216    write (std_out,*) 'vel:'
217    do kk=1,ab_mover%natom
218      write (std_out,*) vel(:,kk)
219    end do
220    write (std_out,*) 'strten:'
221    write (std_out,*) strten(1:3),ch10,strten(4:6)
222    write (std_out,*) 'etotal:'
223    write (std_out,*) etotal
224  end if
225 
226 !Save initial values
227  acell0(:)=acell(:)
228  rprimd0(:,:)=rprimd(:,:)
229  ucvol0=ucvol
230 
231 !Get rid of mean force on whole unit cell, but only if no
232 !generalized constraints are in effect
233 !  call fcart2fred(hist%fcart(:,:,hist%ihist),fred_corrected,rprimd,ab_mover%natom)
234 !  if(ab_mover%nconeq==0)then
235 !    amass_tot=sum(ab_mover%amass(:))
236 !    do ii=1,3
237 !      if (ii/=3.or.ab_mover%jellslab==0) then
238 !        favg=sum(fred_corrected(ii,:))/dble(ab_mover%natom)
239 !        fred_corrected(ii,:)=fred_corrected(ii,:)-favg*ab_mover%amass(:)/amass_tot
240 !      end if
241 !    end do
242 !  end if
243 
244 !write(std_out,*) 'isothermal 03'
245 !##########################################################
246 !### 05. Seconde half velocity step
247 
248  if (itime>1) then
249 
250 !  Next Half velocity step
251    do idim=1,3
252      fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:)
253    end do
254    vel(:,:)=vel_nexthalf(:,:)+ab_mover%dtion/two*fcart_m(:,:)
255 
256    if (ab_mover%optcell==0) then
257 !    Update Thermostat variables and velocity
258      call isotemp(ab_mover%amass,ab_mover%dtion,ekin,ab_mover%iatfix,&
259 &     ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,vel)
260    else if (ab_mover%optcell==1) then
261 !    Update Thermostat variables and velocity
262      call isopress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,&
263 &     ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,&
264 &     strten,ab_mover%strtarget,ucvol,vel,vlogv)
265    else if (ab_mover%optcell==2) then
266 !    Next half step for extended variables
267      call isostress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,&
268 &     ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,&
269 &     ab_mover%qmass,strten,ab_mover%strtarget,ucvol,vel)
270    end if
271 
272    if(itime==2) massvol=ekin+etotal
273 
274    if (ab_mover%optcell==2) then
275 !    Evolution of cell and volume
276      acell_next(:)=acell(:)
277      ucvol_next=ucvol
278    end if
279 
280    call mkrdim(acell,rprim,rprimd)
281    call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
282 
283    if(zDEBUG)then
284      write(std_out,*) 'Second half velocity step'
285      write(std_out,*) 'Cell parameters:'
286      write(std_out,*) 'rprimd:'
287      do kk=1,3
288        write(std_out,*) rprimd(:,kk)
289      end do
290      write(std_out,*) 'rprim:'
291      do kk=1,3
292        write(std_out,*) rprim(:,kk)
293      end do
294      write(std_out,*) 'acell:'
295      write(std_out,*) acell(:)
296      write(std_out,*) 'Conserved energy:',(ekin+etotal)-massvol,ekin,etotal
297      write(std_out,*) 'Volume of unitqry cell (ucvol):',ucvol
298    end if
299 
300  end if ! if (itime>1)
301 
302 !write(std_out,*) 'isothermal 04'
303 !##########################################################
304 !### 03. Compute the next values
305 
306 !The temperature is linear between initial and final values
307 !It is here converted from Kelvin to Hartree (kb_HaK)
308  ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK
309 
310  if(zDEBUG)then
311    write(std_out,*) 'Temperature in Kelvin (ktemp):',ktemp
312    write(std_out,*) 'Initial temp (mdtemp(1)):',ab_mover%mdtemp(1)
313    write(std_out,*) 'Final temp (mdtemp(2)):',ab_mover%mdtemp(2)
314    write(std_out,*) 'Delay for atom permutation (delayperm)',ab_mover%delayperm
315    write(std_out,*) 'nnos:', ab_mover%nnos
316    write(std_out,*) 'qmass', ab_mover%qmass(:)
317    write(std_out,*) 'bmass',ab_mover%bmass
318  end if
319 
320  if(itime==1) then
321    mttk_vars%glogs(:)=zero; mttk_vars%vlogs(:)=zero; mttk_vars%xlogs(:)=zero
322    mttk_vars%vboxg(:,:)=zero
323    vlogv=zero
324 !  v2gauss is twice the kinetic energy
325    v2gauss=0.0_dp
326    do iatom=1,ab_mover%natom
327      do idim=1,3
328        v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
329      end do
330    end do
331 
332 !  If there is no kinetic energy
333    if (v2gauss<=v2tol.and.itime==1) then
334 !    Maxwell-Boltzman distribution
335      v2gauss=zero
336      vtest=zero
337      do iatom=1,ab_mover%natom
338        do idim=1,3
339          vel(idim,iatom)=sqrt(kb_HaK*ab_mover%mdtemp(1)/ab_mover%amass(iatom))*cos(two_pi*uniformrandom(idum))
340          vel(idim,iatom)=vel(idim,iatom)*sqrt(-2._dp*log(uniformrandom(idum)))
341        end do
342      end do
343 
344 !    Get rid of center-of-mass velocity
345      s1=sum(ab_mover%amass(:))
346      do idim=1,3
347        s2=sum(ab_mover%amass(:)*vel(idim,:))
348        vel(idim,:)=vel(idim,:)-s2/s1
349      end do
350 
351 !    Recompute v2gauss
352      do iatom=1,ab_mover%natom
353        do idim=1,3
354          v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
355          vtest=vtest+vel(idim,iatom)/(3._dp*ab_mover%natom)
356        end do
357      end do
358 
359 !    Now rescale the velocities to give the exact temperature
360      rescale_vel=sqrt(3._dp*ab_mover%natom*kb_HaK*ab_mover%mdtemp(1)/v2gauss)
361      vel(:,:)=vel(:,:)*rescale_vel
362 
363 !    Recompute v2gauss with the rescaled velocities
364      v2gauss=zero
365      do iatom=1,ab_mover%natom
366        do idim=1,3
367          v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
368        end do
369      end do
370 
371 !    Compute the variance and print
372      sigma2=(v2gauss/(3._dp*ab_mover%natom)-ab_mover%amass(1)*vtest**2)/kb_HaK
373 
374      if (zDEBUG)then
375        write(message, '(a)' )&
376 &       ' Rescaling or initializing velocities to initial temperature'
377        call wrtout(std_out,message,'COLL')
378        write(message, '(a,d12.5,a,D12.5)' )&
379 &       ' --- Scaling factor :',rescale_vel,' Asked T (K) ',ab_mover%mdtemp(1)
380        call wrtout(std_out,message,'COLL')
381        write(message, '(a,d12.5,a,D12.5)' )&
382 &       ' --- Effective temperature',v2gauss/(3*ab_mover%natom*kb_HaK),' From variance', sigma2
383        call wrtout(std_out,message,'COLL')
384      end if
385 
386    end if !(v2gauss<=v2tol.and.itime==1)
387  end if !(itime==1)
388 
389 !XG070613 : Do not take away the following line , seems needed for the pathscale compiler
390 
391  if (zDEBUG) write(std_out,*) 'vboxg',mttk_vars%vboxg(:,:)
392 
393 
394 !write(std_out,*) 'isothermal 05'
395 !##########################################################
396 !### 03. First half velocity step
397 
398 !write(std_out,*) 'FIRST HALF VELOCITY STEP',ucvol
399 !write(std_out,*) 'OPTCELL option selected:',ab_mover%optcell
400 !write(std_out,*) 'RPRIMD'
401 !do kk=1,3
402 !write(std_out,*) rprimd(:,kk)
403 !end do
404 !write(std_out,*) 'RPRIM'
405 !do kk=1,3
406 !write(std_out,*) rprim(:,kk)
407 !end do
408 !write(std_out,*) 'ACELL'
409 !write(std_out,*) acell(:)
410 
411 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
412 !%%% BEGIN sub case optcell=0 Isothermal Ensemble
413 !%%%
414  if(ab_mover%optcell==0) then
415 !  There is no evolution of cell
416    acell_next(:)=acell(:)
417    ucvol_next=ucvol
418    rprim_next(:,:)=rprim(:,:)
419    rprimd_next(:,:)=rprimd(:,:)
420 !  Update Thermostat variables and scale velocitie
421    call isotemp(ab_mover%amass,ab_mover%dtion,ekin,ab_mover%iatfix,&
422 &   ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,vel)
423 
424 !  Half velocity step
425    do idim=1,3
426      fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:)
427    end do
428    vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:)
429 !  New positions
430 !  Convert input xred (reduced coordinates) to xcart (cartesian)
431    call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
432    xcart_next(:,:)=xcart(:,:)+vel_nexthalf(:,:)*ab_mover%dtion
433 !  Convert back to xred (reduced coordinates)
434    call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
435 !  %%%
436 !  %%% END sub case optcell=0 Isothermal Ensemble
437 !  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
438 
439 !  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
440 !  %%% BEGIN sub case optcell=1 Isothermal-Isenthalpic
441 !  %%%       Ensemble (homogeneous cell deformation)
442 !  %%%
443  else if (ab_mover%optcell==1) then
444 !  Only homogeneous evolution of cell
445 !  Evolution of cell we keep rprim constant
446    rprim_next(:,:)=rprim(:,:)
447 !  Update Thermostat variables and velocity
448    call isopress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,&
449 &   ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,&
450 &   strten,ab_mover%strtarget,ucvol,vel,vlogv)
451 
452 !  Half velocity step
453    do idim=1,3
454      fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:)
455    end do
456    vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:)
457 !  New positions
458    mttk_aloc=exp(ab_mover%dtion/two*vlogv)
459    mttk_aloc2=(vlogv*ab_mover%dtion/two)**2
460    polysh=(((esh8*mttk_aloc2+esh6)*mttk_aloc2+esh4)*mttk_aloc2+esh2)*mttk_aloc2+one
461    mttk_bloc=mttk_aloc*polysh*ab_mover%dtion
462 !  Convert input xred (reduced coordinates) to xcart (cartesian)
463    call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
464    xcart_next(:,:)=xcart(:,:)*mttk_aloc**2+vel_nexthalf(:,:)*mttk_bloc
465 !  Update the volume and related quantities
466    acell_next(:)=acell(:)*exp(ab_mover%dtion*vlogv)
467 !  ucvol=ucvol*exp(ab_mover%dtion*vlogv)
468    call mkrdim(acell_next,rprim,rprimd_next)
469    call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol_next)
470 !  Convert back to xred (reduced coordinates)
471    call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next)
472 !  Computation of the forces for the new positions
473 !  Compute LDA forces (big loop)
474 
475 !  COMMENTED
476 !  This should be in mover.F90
477 
478 !  !      If metric has changed since the initialization, update the Ylm's
479 !  if (ab_mover%optcell/=0.and.psps%useylm==1.and.itime>1)then
480 !  call status(0,dtfil%filstat,iexit,level,'call initylmg ')
481 !  option=0;if (ab_mover%iscf>0) option=1
482 !  call initylmg(gprimd,kg,ab_mover%kptns,ab_mover%mkmem,mpi_enreg,psps%mpsang,ab_mover%mpw,ab_mover%nband,ab_mover%nkpt,&
483 !  &         npwarr,ab_mover%nsppol,option,rprimd_next,ylm,ylmgr)
484 !  end if
485 
486 
487 !  %%%
488 !  %%% END sub case optcell=1 Isothermal-Isenthalpic
489 !  %%%     Ensemble (homogeneous cell deformation)
490 !  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
491 
492 !  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
493 !  %%% BEGIN sub case optcell=2 Isothermal-Isenthalpic
494 !  %%%       Ensemble (full cell deformation)
495 !  %%%
496  else if (ab_mover%optcell==2) then
497    acell_next=acell
498 !  Fisrt half step for extended variables
499    call isostress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,&
500 &   ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,&
501 &   ab_mover%qmass,strten,ab_mover%strtarget,ucvol,vel)
502 !  Half velocity step
503    do idim=1,3
504      fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:)
505    end do
506    vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:)
507 !  Convert input xred (reduced coordinates) to xcart (cartesian)
508    call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
509 !  New positions
510    mttk_vt(:,:)=mttk_vars%vboxg(:,:)
511    call dsyev('V','U',3,mttk_vt,3,mttk_veig,work,lwork,ierr)
512    mttk_tv(:,:)=transpose(mttk_vt)
513    mttk_alc(:)=exp(ab_mover%dtion/two*mttk_veig(:))
514    mttk_alc2(:)=(mttk_veig(:)*ab_mover%dtion/two)**2
515    mttk_psh(:)=(((esh8*mttk_alc2(:)+esh6)*mttk_alc2(:)+esh4)*mttk_alc2(:)+esh2)*mttk_alc2(:)+one
516    mttk_blc(:)=mttk_alc(:)*mttk_psh(:)*ab_mover%dtion
517 !  Update the positions
518    do iatom=1,ab_mover%natom
519      mttk_uu(:)=matmul(mttk_tv,xcart(:,iatom))
520      mttk_uv(:)=matmul(mttk_tv,vel_nexthalf(:,iatom))
521      mttk_uu(:)=mttk_uu(:)*mttk_alc(:)**2+mttk_uv(:)*mttk_blc(:)
522      xcart_next(:,iatom)=matmul(mttk_vt,mttk_uu)
523    end do
524 !  Update the box (rprimd and rprim)
525    mttk_ubox(:,:)=matmul(mttk_tv,rprimd)
526    do idim=1,3
527      mttk_ubox(:,idim)=mttk_ubox(:,idim)*mttk_alc(:)**2
528    end do
529    rprimd_next(:,:)=matmul(mttk_vt,mttk_ubox)
530    do idim=1,3
531      rprim_next(idim,:)=rprimd_next(idim,:)/acell(:)
532    end do
533 !  Update the volume
534    call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol)
535 !  Convert back to xred (reduced coordinates)
536    call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next)
537 !  Computation of the forces for the new positions
538 
539 !  COMMENTED
540 !  This should be in mover.F90
541 
542 !  !      If metric has changed since the initialization, update the Ylm's
543 !  if (ab_mover%optcell/=0.and.psps%useylm==1.and.itime>1)then
544 !  call status(0,dtfil%filstat,iexit,level,'call initylmg ')
545 !  option=0;if (ab_mover%iscf>0) option=1
546 !  call initylmg(gprimd,kg,ab_mover%kptns,ab_mover%mkmem,mpi_enreg,psps%mpsang,ab_mover%mpw,ab_mover%nband,ab_mover%nkpt,&
547 !  &         npwarr,ab_mover%nsppol,option,rprimd_next,ylm,ylmgr)
548 !  end if
549 
550 !  %%%
551 !  %%% END sub case optcell=2 Isothermal-Isenthalpic
552 !  %%%     Ensemble (full cell deformation)
553 !  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
554  else
555    write(message, '(a,i12,a,a)' )&
556 &   '  Disallowed value for optcell=',ab_mover%optcell,ch10,&
557 &   '  Allowed values with ionmov==13 : 0 to 2.'
558    MSG_BUG(message)
559  end if
560 
561 !write(std_out,*) 'OLD PARAMETERS'
562 !write(std_out,*) 'RPRIMD'
563 !do kk=1,3
564 !write(std_out,*) rprimd(:,kk)
565 !end do
566 !write(std_out,*) 'RPRIM'
567 !do kk=1,3
568 !write(std_out,*) rprim(:,kk)
569 !end do
570 !write(std_out,*) 'ACELL'
571 !write(std_out,*) acell(:)
572 
573 !write(std_out,*) 'NEXT PARAMETERS'
574 !write(std_out,*) 'RPRIMD'
575 !do kk=1,3
576 !write(std_out,*) rprimd_next(:,kk)
577 !end do
578 !write(std_out,*) 'RPRIM'
579 !do kk=1,3
580 !write(std_out,*) rprim_next(:,kk)
581 !end do
582 !write(std_out,*) 'ACELL'
583 !write(std_out,*) acell_next(:)
584 
585 
586 !Those are the values store into the history
587  rprim=rprim_next
588  rprimd=rprimd_next
589  xred=xred_next
590  xcart=xcart_next
591  acell=acell_next
592 
593 !write(std_out,*) 'isothermal 06'
594 !##########################################################
595  !### 06. Update the history with the prediction
596 
597 !Increase indexes
598  hist%ihist = abihist_findIndex(hist,+1)
599 
600 !Fill the history with the variables
601 !xred, acell, rprimd, vel
602  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
603  hist%vel(:,:,hist%ihist)=vel(:,:)
604  hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
605 
606  if(zDEBUG)then
607    write (std_out,*) 'fcart:'
608    do kk=1,ab_mover%natom
609      write (std_out,*) fcart(:,kk)
610    end do
611    write (std_out,*) 'vel:'
612    do kk=1,ab_mover%natom
613      write (std_out,*) vel(:,kk)
614    end do
615    write (std_out,*) 'strten:'
616    write (std_out,*) strten(1:3),ch10,strten(4:6)
617    write (std_out,*) 'etotal:'
618    write (std_out,*) etotal
619  end if
620 
621 end subroutine pred_isothermal