diff --git a/Src/config/make.inc b/Src/config/make.inc index a38eb60..1509b39 100644 --- a/Src/config/make.inc +++ b/Src/config/make.inc @@ -20,12 +20,14 @@ CPU=Intel # Suggestions for compilers. Feel free to override the suggestions. ifeq ($(CPU),Intel) - F77=ifort + F77=ifx + # F77=ifort # F77=gfortran # F77=pgfortran else ifeq ($(CPU),AMD) F77=gfortran + # F77=ifx # F77=ifort # F77=pgfortran else @@ -68,9 +70,17 @@ PRE2008=false -# For Intel ifort compiler: +# For Intel ifort and ifx compiler: +IFORTX=false ifeq ($(F77),ifort) + IFORTX=true +endif +ifeq ($(F77),ifx) + IFORTX=true +endif + +ifeq ($(IFORTX),true) # profiling option (-g,empty). Use it for tunning the code. @@ -116,7 +126,7 @@ ifeq ($(F77),ifort) FCOMP=-qopenmp # FCOMP=-openmp else - FCOPTS=-check all -warn all,nodec,interfaces,noexternal -gen-interfaces -traceback -fpe0 -fp-stack-check -g + FCOPTS=-check all -warn all,nodec,interfaces,noexternal -gen-interfaces -traceback -fpe0 -fp-stack-check -g -O0 FCOMP= endif diff --git a/Src/cpw.f90 b/Src/cpw.f90 index 58df9b2..a4076bc 100644 --- a/Src/cpw.f90 +++ b/Src/cpw.f90 @@ -124,7 +124,7 @@ program cpw2000 ! Driver program version - vdriv = '5.10' + vdriv = '5.11' ! timing diff --git a/Src/cpw/cpw_force.f90 b/Src/cpw/cpw_force.f90 index 7370c09..158190a 100644 --- a/Src/cpw/cpw_force.f90 +++ b/Src/cpw/cpw_force.f90 @@ -11,119 +11,143 @@ ! https://github.com/jlm785/cpw2000 ! !------------------------------------------------------------! -!> manages the calculation of forces and stresses. -!> it is an interface to force_stress_kb with some extra functions - - subroutine cpw_force(iprglob,strxc, ealpha, deltentpy, errfrc, & - & errstr, & - & dims_, crys_, total_, ewald_, flags_, spaceg_, recip_, pseudo_, & - & vcomp_, chdens_, hamallk_, psiallk_, kpoint_, vcsdyn_) - -! Written November 2019. JLM -! Modified upstream in January 2010. vff_add_keating. JLM -! Modified, error after keating, 18 February 2020. JLM -! copyright J.L.Martins, INESC-MN. - -! version 4.95 +!> Manages the calculation of forces and stresses. +!> it is an interface to force_stress_kb with the extra functions +!> of adding the Keating corrections and printing the result. +!> +!> \author Jose Luis Martins +!> \version 5.10 +!> \date 20 October 93, 21 February 2024. +!> \copyright GNU Public License v2 + +subroutine cpw_force(iprglob,strxc, ealpha, deltentpy, errfrc, & + errstr, & + dims_, crys_, total_, ewald_, flags_, spaceg_, recip_, pseudo_, & + vcomp_, chdens_, hamallk_, psiallk_, kpoint_, vcsdyn_) + +! Written November 2019. JLM +! Modified upstream in January 2010. vff_add_keating. JLM +! Modified, error after keating, 18 February 2020. JLM +! Modified, indentation, another printing choices, 21 February 2024. JLM + + + + use cpw_variables + + implicit none + + type(dims_t) :: dims_ !< array dimensions + type(crys_t) :: crys_ !< crystal structure + type(enfrst_t) :: total_ !< Total energy force stress + type(enfrst_t) :: ewald_ !< Ewald energy force stress + type(flags_t) :: flags_ !< computational flags + type(spaceg_t) :: spaceg_ !< space group information + type(recip_t) :: recip_ !< reciprocal space information + type(pseudo_t) :: pseudo_ !< pseudo-potential (Kleinman-Bylander) + type(vcomp_t) :: vcomp_ !< Componemts of local potential + type(chdens_t) :: chdens_ !< charge densities + type(hamallk_t) :: hamallk_ !< hamiltonian size and indexation for all k-points + type(psiallk_t) :: psiallk_ !< psi for all k-points + type(kpoint_t) :: kpoint_ !< k-point data + type(vcsdyn_t) :: vcsdyn_ !< variational cell shape molecular dynamics variables + + integer, intent(in) :: iprglob !< controls the amount of printing by subroutines + real(REAL64), intent(inout) :: strxc(3,3) !< contribution of xc to the stress tensor (contravariant,Hartree) + real(REAL64), intent(in) :: ealpha !< G=0 contribution to the total energy (Hartree) + + real(REAL64), intent(out) :: deltentpy !< enthalpy difference from last iteration + real(REAL64), intent(out) :: errfrc !< maximum error in force (cartesian coordiantes) Hartree/Bohr + real(REAL64), intent(out) :: errstr !< maximum error in stress + + integer :: ipr + integer :: iotape + + integer :: minrat ! if =1 minimize with respect to atomic positions + integer :: minstr ! if =1 minimize with respect to all adot variables, if =2 minimize with respect to adot(3,3) + + call force_stress_kb(total_%force, total_%stress, total_%energy, & + ewald_%force, ewald_%stress, strxc, ealpha, flags_%flgpsd, & + crys_%ntype, crys_%natom, crys_%nameat, crys_%rat, crys_%adot, & + spaceg_%ntrans, spaceg_%mtrx, spaceg_%tnp, & + recip_%ng, recip_%kgv, recip_%phase, recip_%conj, recip_%ns, & + recip_%mstar, recip_%ek, & + pseudo_%nq, pseudo_%delq, pseudo_%vkb, pseudo_%nkb, & + vcomp_%vion, vcomp_%vhar, vcomp_%vxc, chdens_%den, & + hamallk_%mtxd_allk, hamallk_%isort_allk, & + psiallk_%psi_allk, psiallk_%occ_allk, & + pseudo_%vql, pseudo_%dnc, pseudo_%dvql, pseudo_%ddc, & + kpoint_%nrk, kpoint_%nband, kpoint_%rk, kpoint_%wgk, & + dims_%mxdtyp, dims_%mxdatm, dims_%mxdlqp, dims_%mxddim, & + dims_%mxdbnd, dims_%mxdgve, dims_%mxdnst, dims_%mxdnrk) + + +!-----------------------clr---------------------------------------- + + + ipr = 0 + if(iprglob > 1) ipr = 1 + + call print_energy(ipr, 'Total', & + total_%energy, total_%force, total_%stress, & + crys_%adot, crys_%ntype, crys_%natom, crys_%nameat, & + dims_%mxdtyp, dims_%mxdatm) + +! Adds keating correction + + if(flags_%flgkeat == 'KEATNG') then + + write(6,*) + write(6,*) ' Using Keating Force Field Corrections' + write(6,*) + + ipr = 0 + if(iprglob > 2) ipr = 1 + iotape = 6 + + call vff_add_keating(ipr, iotape, & + crys_%ntype, crys_%natom, crys_%nameat, & + crys_%rat, crys_%adot, & + total_%force, total_%energy, total_%stress, & + dims_%mxdtyp, dims_%mxdatm) + + ipr = 1 + call print_energy(ipr, 'Total+Keating', & + total_%energy, total_%force, total_%stress, & + crys_%adot, crys_%ntype, crys_%natom, crys_%nameat, & + dims_%mxdtyp, dims_%mxdatm) + + else + ipr = 0 + if(iprglob < 2) ipr = 1 - use cpw_variables - - implicit none - - type(dims_t) :: dims_ !< array dimensions - type(crys_t) :: crys_ !< crystal structure - type(enfrst_t) :: total_ !< Total energy force stress type(flags_t) :: flags_ !< computational flags - type(enfrst_t) :: ewald_ !< Ewald energy force stress - type(flags_t) :: flags_ !< computational flags - type(spaceg_t) :: spaceg_ !< space group information - type(recip_t) :: recip_ !< reciprocal space information - type(pseudo_t) :: pseudo_ !< pseudo-potential (Kleinman-Bylander) - type(vcomp_t) :: vcomp_ !< Componemts of local potential - type(chdens_t) :: chdens_ !< charge densities - type(hamallk_t) :: hamallk_ !< hamiltonian size and indexation for all k-points - type(psiallk_t) :: psiallk_ !< psi for all k-points - type(kpoint_t) :: kpoint_ !< k-point data - type(vcsdyn_t) :: vcsdyn_ !< variational cell shape molecular dynamics variables - - integer, intent(in) :: iprglob !< controls the amount of printing by subroutines - real(REAL64), intent(inout) :: strxc(3,3) !< contribution of xc to the stress tensor (contravariant,Hartree) - real(REAL64), intent(in) :: ealpha !< G=0 contribution to the total energy (Hartree) - - real(REAL64), intent(out) :: deltentpy !< enthalpy difference from last iteration - real(REAL64), intent(out) :: errfrc !< maximum error in force (cartesian coordiantes) Hartree/Bohr - real(REAL64), intent(out) :: errstr !< maximum error in stress - - integer :: ipr - integer :: iotape - - integer :: minrat ! if =1 minimize with respect to atomic positions - integer :: minstr ! if =1 minimize with respect to all adot variables, if =2 minimize with respect to adot(3,3) - - call force_stress_kb(total_%force, total_%stress, total_%energy, & - & ewald_%force, ewald_%stress, strxc, ealpha, flags_%flgpsd, & - & crys_%ntype, crys_%natom, crys_%nameat, crys_%rat, crys_%adot, & - & spaceg_%ntrans, spaceg_%mtrx, spaceg_%tnp, & - & recip_%ng, recip_%kgv, recip_%phase, recip_%conj, recip_%ns, & - & recip_%mstar, recip_%ek, & - & pseudo_%nq, pseudo_%delq, pseudo_%vkb, pseudo_%nkb, & - & vcomp_%vion, vcomp_%vhar, vcomp_%vxc, chdens_%den, & - & hamallk_%mtxd_allk, hamallk_%isort_allk, & - & psiallk_%psi_allk, psiallk_%occ_allk, & - & pseudo_%vql, pseudo_%dnc, pseudo_%dvql, pseudo_%ddc, & - & kpoint_%nrk, kpoint_%nband, kpoint_%rk, kpoint_%wgk, & - & dims_%mxdtyp, dims_%mxdatm, dims_%mxdlqp, dims_%mxddim, & - & dims_%mxdbnd, dims_%mxdgve, dims_%mxdnst, dims_%mxdnrk) - - -!----------------------------clr---------------------------------------- - -! Adds keating correction - - if(flags_%flgkeat == 'KEATNG') then - - write(6,*) 'Using Keating ForceField Corrections' - - ipr = 0 - if(iprglob > 2) ipr = 1 - iotape = 6 - - call vff_add_keating(ipr, iotape, & - & crys_%ntype, crys_%natom, crys_%nameat, & - & crys_%rat, crys_%adot, & - & total_%force, total_%energy, total_%stress, & - & dims_%mxdtyp, dims_%mxdatm) - - ipr = 0 - if(iprglob > 1) ipr = 1 - - call print_energy(ipr, 'Total+Keating', & - & total_%energy, total_%force, total_%stress, & - & crys_%adot, crys_%ntype, crys_%natom, crys_%nameat, & - & dims_%mxdtyp, dims_%mxdatm) - - endif - - if(flags_%flgcal == 'LBFSYM' .or. flags_%flgcal == 'VCSLBF' .or. & - & flags_%flgcal == 'EPILBF') then - - minrat = 1 - if(flags_%flgcal == 'LBFSYM') then - minstr = 0 - elseif(flags_%flgcal == 'VCSLBF') then - minstr = 1 - elseif(flags_%flgcal == 'EPILBF') then - minstr = 2 - endif - - call force_stress_error(total_%energy, total_%force, & - & total_%stress, vcsdyn_%press, vcsdyn_%strext, & - & deltentpy, errfrc, errstr, minrat, minstr, & - & crys_%adot, crys_%ntype, crys_%natom, & - & dims_%mxdtyp, dims_%mxdatm) - - endif + call print_energy(ipr, 'Total', & + total_%energy, total_%force, total_%stress, & + crys_%adot, crys_%ntype, crys_%natom, crys_%nameat, & + dims_%mxdtyp, dims_%mxdatm) - return + endif - end subroutine cpw_force + if(flags_%flgcal == 'LBFSYM' .or. flags_%flgcal == 'VCSLBF' .or. & + flags_%flgcal == 'EPILBF') then + + minrat = 1 + if(flags_%flgcal == 'LBFSYM') then + minstr = 0 + elseif(flags_%flgcal == 'VCSLBF') then + minstr = 1 + elseif(flags_%flgcal == 'EPILBF') then + minstr = 2 + endif + + call force_stress_error(total_%energy, total_%force, & + total_%stress, vcsdyn_%press, vcsdyn_%strext, & + deltentpy, errfrc, errstr, minrat, minstr, & + crys_%adot, crys_%ntype, crys_%natom, & + dims_%mxdtyp, dims_%mxdatm) + + endif + + return + +end subroutine cpw_force diff --git a/Src/cpw_lib/force_stress_kb.f90 b/Src/cpw_lib/force_stress_kb.f90 index 129c50d..7f95fa2 100644 --- a/Src/cpw_lib/force_stress_kb.f90 +++ b/Src/cpw_lib/force_stress_kb.f90 @@ -11,324 +11,321 @@ ! https://github.com/jlm785/cpw2000 ! !------------------------------------------------------------! -!> Calculates the force and stress - - subroutine force_stress_kb(force,stress,energy, & - & forcew,strew,strxc,ealpha,flgpsd, & - & ntype,natom,nameat,rat,adot, & - & ntrans,mtrx,tnp, & - & ng,kgv,phase,conj,ns,mstar,ek, & - & nqnl,delqnl,VKB,nkb, & - & vion,vhar,vxc,den, & - & mtxd_allk,isort_allk,psi_allk,occ_allk, & - & vql,dnc,dvql,ddc, & - & nrk,nband,rk,wgk, & - & mxdtyp,mxdatm,mxdlqp,mxddim,mxdbnd,mxdgve,mxdnst,mxdnrk) - -! Calculates the force and stress - -! Version 4.0. 20 october 93. jlm -! Version 4.40 24/07/2002. jlm -! Modified to symmetrize strxc (gga) 28/9/2004, version 4.42 jlm -! Modified f90, January 2017. JLM -! Modified, documentation, January 2020. JLM -! copyright inesc-mn/Jose Luis Martins - -! version 4.94 +!> Calculates the force and stress +!> +!> \author Jose Luis Martins +!> \version 5.10 +!> \date 20 October 93, 21 February 2024. +!> \copyright GNU Public License v2 + + subroutine force_stress_kb(force, stress, energy, & + forcew, strew, strxc, ealpha, flgpsd, & + ntype, natom, nameat, rat, adot, & + ntrans, mtrx, tnp, & + ng, kgv, phase, conj, ns, mstar, ek, & + nqnl, delqnl, vkb, nkb, & + vion, vhar, vxc, den, & + mtxd_allk, isort_allk, psi_allk, occ_allk, & + vql, dnc, dvql, ddc, & + nrk, nband, rk, wgk, & + mxdtyp, mxdatm, mxdlqp, mxddim, mxdbnd, mxdgve, mxdnst, mxdnrk) + +! Version 4.0. 20 october 93. jlm +! Version 4.40 24/07/2002. jlm +! Modified to symmetrize strxc (gga) 28/9/2004, version 4.42 jlm +! Modified f90, January 2017. JLM +! Modified, documentation, January 2020. JLM +! Modified, indentation, remove print, 21 February 2024. JLM + - implicit none + implicit none - integer, parameter :: REAL64 = selected_real_kind(12) + integer, parameter :: REAL64 = selected_real_kind(12) -! input +! input - integer, intent(in) :: mxdtyp !< array dimension of types of atoms - integer, intent(in) :: mxdatm !< array dimension of number of atoms of a given type - integer, intent(in) :: mxdlqp !< array dimension for local potential - integer, intent(in) :: mxddim !< array dimension of plane-waves - integer, intent(in) :: mxdbnd !< array dimension for number of bands - integer, intent(in) :: mxdgve !< array dimension for g-space vectors - integer, intent(in) :: mxdnst !< array dimension for g-space stars - integer, intent(in) :: mxdnrk !< array dimension of k-points + integer, intent(in) :: mxdtyp !< array dimension of types of atoms + integer, intent(in) :: mxdatm !< array dimension of number of atoms of a given type + integer, intent(in) :: mxdlqp !< array dimension for local potential + integer, intent(in) :: mxddim !< array dimension of plane-waves + integer, intent(in) :: mxdbnd !< array dimension for number of bands + integer, intent(in) :: mxdgve !< array dimension for g-space vectors + integer, intent(in) :: mxdnst !< array dimension for g-space stars + integer, intent(in) :: mxdnrk !< array dimension of k-points - real(REAL64), intent(in) :: energy !< Total electronic energy (Hartree) - real(REAL64), intent(in) :: forcew(3,mxdatm,mxdtyp) !< d enerew / d rat, Ewald contribution to force (contravariant components) - real(REAL64), intent(in) :: strew(3,3) !< d enerew / d adot, Ewald contribution to stress tensor (contravariant components) - real(REAL64), intent(in) :: ealpha !< G=0 contribution to the total energy (Hartree) - character(len=6), intent(in) :: flgpsd !< type of pseudopotential + real(REAL64), intent(in) :: energy !< Total electronic energy (Hartree) + real(REAL64), intent(in) :: forcew(3,mxdatm,mxdtyp) !< d enerew / d rat, Ewald contribution to force (contravariant components) + real(REAL64), intent(in) :: strew(3,3) !< d enerew / d adot, Ewald contribution to stress tensor (contravariant components) + real(REAL64), intent(in) :: ealpha !< G=0 contribution to the total energy (Hartree) + character(len=6), intent(in) :: flgpsd !< type of pseudopotential - integer, intent(in) :: ntype !< number of types of atoms - integer, intent(in) :: natom(mxdtyp) !< number of atoms of type i - character(len=2), intent(in) :: nameat(mxdtyp) !< chemical symbol for the type i - real(REAL64), intent(in) :: rat(3,mxdatm,mxdtyp) !< lattice coordinates of atom j of type i - real(REAL64), intent(in) :: adot(3,3) !< metric in real space + integer, intent(in) :: ntype !< number of types of atoms + integer, intent(in) :: natom(mxdtyp) !< number of atoms of type i + character(len=2), intent(in) :: nameat(mxdtyp) !< chemical symbol for the type i + real(REAL64), intent(in) :: rat(3,mxdatm,mxdtyp) !< lattice coordinates of atom j of type i + real(REAL64), intent(in) :: adot(3,3) !< metric in real space - integer, intent(in) :: ntrans !< number of symmetry operations in the factor group - integer, intent(in) :: mtrx(3,3,48) !< rotation matrix (in reciprocal lattice coordinates) for the k-th symmetry operation of the factor group - real(REAL64), intent(in) :: tnp(3,48) !< 2*pi* i-th component (in lattice coordinates) of the fractional translation vector associated with the k-th symmetry operation of the factor group - - integer, intent(in) :: ng !< total number of g-vectors with length less than gmax - integer, intent(in) :: kgv(3,mxdgve) !< i-th component (reciprocal lattice coordinates) of the n-th g-vector ordered by stars of increasing length - complex(REAL64), intent(in) :: phase(mxdgve) !< real part of the phase factor of G-vector n - real(REAL64), intent(in) :: conj(mxdgve) !< is -1 if one must take the complex conjugate of x*phase - integer, intent(in) :: ns !< number os stars with length less than gmax - integer, intent(in) :: mstar(mxdnst) !< number of G-vectors in the j-th star - real(REAL64), intent(in) :: ek(mxdnst) !< kinetic energy (hartree) of g-vectors in star j - - integer, intent(in) :: nqnl(mxdtyp) !< number of points for the non-local pseudopotential interpolation - real(REAL64), intent(in) :: delqnl(mxdtyp) !< step used in the interpolation - real(REAL64), intent(in) :: vkb(-2:mxdlqp,0:3,-1:1,mxdtyp) !< (1/q**l) * KB nonlocal pseudo. for atom k, ang. mom. l. normalized to vcell, hartree - integer, intent(in) :: nkb(0:3,-1:1,mxdtyp) !< KB pseudo. normalization for atom k, ang. mom. l - - complex(REAL64), intent(in) :: vion(mxdnst) !< ionic potential for the prototype G-vector in star j - complex(REAL64), intent(in) :: den(mxdnst) !< total charge density for the prototype G-vector - complex(REAL64), intent(in) :: vhar(mxdnst) !< Hartree potential for the prototype G-vector - complex(REAL64), intent(in) :: vxc(mxdnst) !< exchange+correlation potential for the prototype G-vector - - integer, intent(in) :: mtxd_allk(mxdnrk) !< dimension of the hamiltonian for k-point n - integer, intent(in) :: isort_allk(mxddim,mxdnrk) !< G-vector associated with k+G vector i of hamiltonian for k-point n - complex(REAL64), intent(in) :: psi_allk(mxddim,mxdbnd,mxdnrk) !< eigenvectors for all k-points - real(REAL64), intent(in) :: occ_allk(mxdnrk*mxdbnd) !< fractional ocupation of level j, for all the k-points - - real(REAL64), intent(in) :: vql(mxdtyp,mxdnst) !< local pseudopotential for atom type i and prototype g-vector in star j real*8 floc(3,mxdatm,mxdtyp) - real(REAL64), intent(in) :: dnc(mxdtyp,mxdnst) !< core charge for atom type i and prototype g-vector in star j - complex(REAL64), intent(in) :: dvql(mxdnst) !< derivative of the local pseudopotential for the prototype g-vector in star j - complex(REAL64), intent(in) :: ddc(mxdnst) !< derivative of the core charge for the prototype g-vector in star j - - integer, intent(in) :: nrk !< number of k-points for integration in the irreducible wedge of the brillouin zone - integer, intent(in) :: nband(mxdnrk) !< number of bands for each k-points - real(REAL64), intent(in) :: rk(3,mxdnrk) !< component in lattice coordinates of the k-point in the mesh - real(REAL64), intent(in) :: wgk(mxdnrk) !< weight in the integration of k-point - -! input and output - - real(REAL64), intent(inout) :: strxc(3,3) !< contribution of xc to the stress tensor (contravariant,Hartree). It is only symmetrized here. - -! output - - real(REAL64), intent(out) :: force(3,mxdatm,mxdtyp) !< d energy / d rat, force on the n-th atom of type i (contravarian components, hartree/bohr) - real(REAL64), intent(out) :: stress(3,3) !< d energy / d adot, stress tensor (contravariant) - -! local variables - - real(REAL64) :: vcell, bdot(3,3), adotm1(3,3) - real(REAL64) :: rkpt(3) - real(REAL64) :: sunsym(3,3), ssym(3,3) - real(REAL64) :: strhl(3,3), strnlkb(3,3), strkin(3,3) - integer :: iel, neig, ipr, mtxd - -! allocatable variables - - real(REAL64), allocatable :: occp(:) ! ocupation*weight*spin deg. of eigenvector j - - real(REAL64), allocatable :: floc(:,:,:) - real(REAL64), allocatable :: funsym(:,:,:) - real(REAL64), allocatable :: fsym(:,:,:) - real(REAL64), allocatable :: fnlkb(:,:,:) - -! constants - - real(REAL64), parameter :: ZERO = 0.0_REAL64 - real(REAL64), parameter :: PI = 3.14159265358979323846_REAL64 - -! counters - - integer :: i, j, k, l, irk - - - - allocate(occp(mxdbnd)) - - allocate(floc(3,mxdatm,mxdtyp)) - allocate(funsym(3,mxdatm,mxdtyp)) - allocate(fsym(3,mxdatm,mxdtyp)) - allocate(fnlkb(3,mxdatm,mxdtyp)) - - - if(flgpsd /= 'PSEUKB') then - write(6,'(" STOPPED in force_stress_kb: wrong ", & - & "pseudopotencial ",a6)') flgpsd - - stop - - endif - - call adot_to_bdot(adot,vcell,bdot) - - do i = 1,3 - do j = 1,3 - adotm1(i,j) = bdot(i,j) / (4*PI*PI) - enddo - enddo - -! local contribution to force (covariant) - - call for_str_local_force(floc, & - & ng,kgv,phase,conj,ns,mstar, & - & vxc,den, & - & ntype,natom,rat, & - & vql,dnc, & - & mxdgve,mxdnst,mxdtyp,mxdatm) - -! local contribution to stress (covariant) - - call for_str_local_stress(ealpha,strhl, & - & ng,kgv,ns,mstar,ek, & - & vion,vhar,vxc,den, & - & adot, & - & dvql,ddc, & - & mxdgve,mxdnst) - -! strxc in GGA is not strictly symmetric. -! It is not elegant but each contribution is kept separate - - do i = 1,3 - do j = 1,3 - sunsym(i,j) = ZERO - do k=1,3 - do l=1,3 - sunsym(i,j) = sunsym(i,j) + adot(i,k)*strxc(k,l)*adot(l,j) - enddo - enddo - enddo - enddo - - do i = 1,ntype - do j = 1,natom(i) - do k = 1,3 - funsym(k,j,i) = ZERO - enddo - enddo - enddo - -! symetrizes - - call for_str_sym(funsym,sunsym,fsym,ssym, & - & ntype,natom,rat, & - & ntrans,mtrx,tnp, & - & mxdtyp,mxdatm) - - do i = 1,3 - do j = 1,3 - strxc(i,j) = ZERO - do k = 1,3 - do l = 1,3 - strxc(i,j) = strxc(i,j) + adotm1(i,k)*ssym(k,l)*adotm1(l,j) - enddo - enddo - enddo - enddo - - do i = 1,3 - do j = 1,3 - sunsym(j,i) = ZERO - enddo - enddo - - do i=1,ntype - do j=1,natom(i) - do k=1,3 - funsym(k,j,i) = ZERO - enddo - enddo - enddo - - iel = 0 - do irk = 1,nrk - -! loop over k-points - - rkpt(1) = rk(1,irk) - rkpt(2) = rk(2,irk) - rkpt(3) = rk(3,irk) - - neig = nband(irk) - mtxd = mtxd_allk(irk) - - do j = 1,neig - iel = iel + 1 - occp(j) = 2*wgk(irk)*occ_allk(iel) - enddo -! - call for_str_kinetic_stress(strkin, & - & mtxd,rkpt,neig,occp, & - & isort_allk(:,irk),psi_allk(:,:,irk), & - & kgv, & - & mxdgve,mxddim,mxdbnd) - - do i = 1,3 - do j = 1,3 - sunsym(j,i) = sunsym(j,i) + strkin(j,i) - enddo - enddo - - call for_str_nl_kb(fnlkb,strnlkb, & - & mtxd,rkpt,neig,occp, & - & isort_allk(:,irk),psi_allk(:,:,irk), & - & kgv, & - & nqnl,delqnl,vkb,nkb, & - & ntype,natom,rat,adot, & - & mxdtyp,mxdatm,mxdgve,mxdlqp,mxddim,mxdbnd) - - do i = 1,3 - do j = 1,3 - sunsym(j,i) = sunsym(j,i) + strnlkb(j,i) - enddo - enddo - - do i = 1,ntype - do j=1,natom(i) - do k=1,3 - funsym(k,j,i) = funsym(k,j,i) + fnlkb(k,j,i) - enddo - enddo - enddo - - enddo - - call for_str_sym(funsym,sunsym,fsym,ssym, & - & ntype,natom,rat, & - & ntrans,mtrx,tnp, & - & mxdtyp,mxdatm) - -! sums up contributions - - do i = 1,3 - do j = 1,3 - stress(i,j) = strew(i,j) + strxc(i,j) - do k=1,3 - do l=1,3 - stress(i,j) = stress(i,j) + adotm1(i,k) * & - & (strhl(k,l) + ssym(k,l)) * adotm1(l,j) - enddo - enddo - enddo - enddo - - do i=1,ntype - do j=1,natom(i) - do k=1,3 - force(k,j,i) = forcew(k,j,i) - do l=1,3 - force(k,j,i) = force(k,j,i) + adotm1(k,l) & - & * (fsym(l,j,i) + floc(l,j,i)) - enddo - enddo - enddo - enddo - - ipr = 1 - call print_energy(ipr,'Total',energy,force,stress, & - & adot,ntype,natom,nameat, & - & mxdtyp,mxdatm) - - deallocate(occp) - - deallocate(floc) - deallocate(funsym) - deallocate(fsym) - deallocate(fnlkb) - - return - end subroutine force_stress_kb + integer, intent(in) :: ntrans !< number of symmetry operations in the factor group + integer, intent(in) :: mtrx(3,3,48) !< rotation matrix (in reciprocal lattice coordinates) for the k-th symmetry operation of the factor group + real(REAL64), intent(in) :: tnp(3,48) !< 2*pi* i-th component (in lattice coordinates) of the fractional translation vector associated with the k-th symmetry operation of the factor group + + integer, intent(in) :: ng !< total number of g-vectors with length less than gmax + integer, intent(in) :: kgv(3,mxdgve) !< i-th component (reciprocal lattice coordinates) of the n-th g-vector ordered by stars of increasing length + complex(REAL64), intent(in) :: phase(mxdgve) !< real part of the phase factor of G-vector n + real(REAL64), intent(in) :: conj(mxdgve) !< is -1 if one must take the complex conjugate of x*phase + integer, intent(in) :: ns !< number os stars with length less than gmax + integer, intent(in) :: mstar(mxdnst) !< number of G-vectors in the j-th star + real(REAL64), intent(in) :: ek(mxdnst) !< kinetic energy (hartree) of g-vectors in star j + + integer, intent(in) :: nqnl(mxdtyp) !< number of points for the non-local pseudopotential interpolation + real(REAL64), intent(in) :: delqnl(mxdtyp) !< step used in the interpolation + real(REAL64), intent(in) :: vkb(-2:mxdlqp,0:3,-1:1,mxdtyp) !< (1/q**l) * KB nonlocal pseudo. for atom k, ang. mom. l. normalized to vcell, hartree + integer, intent(in) :: nkb(0:3,-1:1,mxdtyp) !< KB pseudo. normalization for atom k, ang. mom. l + + complex(REAL64), intent(in) :: vion(mxdnst) !< ionic potential for the prototype G-vector in star j + complex(REAL64), intent(in) :: den(mxdnst) !< total charge density for the prototype G-vector + complex(REAL64), intent(in) :: vhar(mxdnst) !< Hartree potential for the prototype G-vector + complex(REAL64), intent(in) :: vxc(mxdnst) !< exchange+correlation potential for the prototype G-vector + + integer, intent(in) :: mtxd_allk(mxdnrk) !< dimension of the hamiltonian for k-point n + integer, intent(in) :: isort_allk(mxddim,mxdnrk) !< G-vector associated with k+G vector i of hamiltonian for k-point n + complex(REAL64), intent(in) :: psi_allk(mxddim,mxdbnd,mxdnrk) !< eigenvectors for all k-points + real(REAL64), intent(in) :: occ_allk(mxdnrk*mxdbnd) !< fractional ocupation of level j, for all the k-points + + real(REAL64), intent(in) :: vql(mxdtyp,mxdnst) !< local pseudopotential for atom type i and prototype g-vector in star j real*8 floc(3,mxdatm,mxdtyp) + real(REAL64), intent(in) :: dnc(mxdtyp,mxdnst) !< core charge for atom type i and prototype g-vector in star j + complex(REAL64), intent(in) :: dvql(mxdnst) !< derivative of the local pseudopotential for the prototype g-vector in star j + complex(REAL64), intent(in) :: ddc(mxdnst) !< derivative of the core charge for the prototype g-vector in star j + + integer, intent(in) :: nrk !< number of k-points for integration in the irreducible wedge of the brillouin zone + integer, intent(in) :: nband(mxdnrk) !< number of bands for each k-points + real(REAL64), intent(in) :: rk(3,mxdnrk) !< component in lattice coordinates of the k-point in the mesh + real(REAL64), intent(in) :: wgk(mxdnrk) !< weight in the integration of k-point + +! input and output + + real(REAL64), intent(inout) :: strxc(3,3) !< contribution of xc to the stress tensor (contravariant,Hartree). It is only symmetrized here. + +! output + + real(REAL64), intent(out) :: force(3,mxdatm,mxdtyp) !< d energy / d rat, force on the n-th atom of type i (contravarian components, hartree/bohr) + real(REAL64), intent(out) :: stress(3,3) !< d energy / d adot, stress tensor (contravariant) + +! local variables + + real(REAL64) :: vcell, bdot(3,3), adotm1(3,3) + real(REAL64) :: rkpt(3) + real(REAL64) :: sunsym(3,3), ssym(3,3) + real(REAL64) :: strhl(3,3), strnlkb(3,3), strkin(3,3) + integer :: iel, neig, ipr, mtxd + +! allocatable variables + + real(REAL64), allocatable :: occp(:) ! ocupation*weight*spin deg. of eigenvector j + + real(REAL64), allocatable :: floc(:,:,:) + real(REAL64), allocatable :: funsym(:,:,:) + real(REAL64), allocatable :: fsym(:,:,:) + real(REAL64), allocatable :: fnlkb(:,:,:) + +! constants + + real(REAL64), parameter :: ZERO = 0.0_REAL64 + real(REAL64), parameter :: PI = 3.14159265358979323846_REAL64 + +! counters + + integer :: i, j, k, l, irk + + + + allocate(occp(mxdbnd)) + + allocate(floc(3,mxdatm,mxdtyp)) + allocate(funsym(3,mxdatm,mxdtyp)) + allocate(fsym(3,mxdatm,mxdtyp)) + allocate(fnlkb(3,mxdatm,mxdtyp)) + + + if(flgpsd /= 'PSEUKB') then + write(6,'(" STOPPED in force_stress_kb: wrong pseudopotencial ",a6)') flgpsd + + stop + + endif + + call adot_to_bdot(adot,vcell,bdot) + + do i = 1,3 + do j = 1,3 + adotm1(i,j) = bdot(i,j) / (4*PI*PI) + enddo + enddo + +! local contribution to force (covariant) + + call for_str_local_force(floc, & + ng, kgv, phase, conj, ns, mstar, & + vxc, den, & + ntype, natom, rat, & + vql, dnc, & + mxdgve, mxdnst, mxdtyp, mxdatm) + +! local contribution to stress (covariant) + + call for_str_local_stress(ealpha, strhl, & + ng, kgv, ns, mstar, ek, & + vion, vhar, vxc, den, & + adot, & + dvql, ddc, & + mxdgve, mxdnst) + +! strxc in GGA is not strictly symmetric. +! It is not elegant but each contribution is kept separate + + do i = 1,3 + do j = 1,3 + sunsym(i,j) = ZERO + do k=1,3 + do l=1,3 + sunsym(i,j) = sunsym(i,j) + adot(i,k)*strxc(k,l)*adot(l,j) + enddo + enddo + enddo + enddo + + do i = 1,ntype + do j = 1,natom(i) + do k = 1,3 + funsym(k,j,i) = ZERO + enddo + enddo + enddo + +! symetrizes + + call for_str_sym(funsym, sunsym, fsym, ssym, & + ntype, natom, rat, & + ntrans, mtrx, tnp, & + mxdtyp ,mxdatm) + + do i = 1,3 + do j = 1,3 + strxc(i,j) = ZERO + do k = 1,3 + do l = 1,3 + strxc(i,j) = strxc(i,j) + adotm1(i,k)*ssym(k,l)*adotm1(l,j) + enddo + enddo + enddo + enddo + + do i = 1,3 + do j = 1,3 + sunsym(j,i) = ZERO + enddo + enddo + + do i = 1,ntype + do j = 1,natom(i) + do k = 1,3 + funsym(k,j,i) = ZERO + enddo + enddo + enddo + + iel = 0 + do irk = 1,nrk + +! loop over k-points + + rkpt(1) = rk(1,irk) + rkpt(2) = rk(2,irk) + rkpt(3) = rk(3,irk) + + neig = nband(irk) + mtxd = mtxd_allk(irk) + + do j = 1,neig + iel = iel + 1 + occp(j) = 2*wgk(irk)*occ_allk(iel) + enddo +! + call for_str_kinetic_stress(strkin, & + mtxd, rkpt, neig, occp, & + isort_allk(:,irk) ,psi_allk(:,:,irk), & + kgv, & + mxdgve, mxddim, mxdbnd) + + do i = 1,3 + do j = 1,3 + sunsym(j,i) = sunsym(j,i) + strkin(j,i) + enddo + enddo + + call for_str_nl_kb(fnlkb, strnlkb, & + mtxd, rkpt, neig, occp, & + isort_allk(:,irk), psi_allk(:,:,irk), & + kgv, & + nqnl, delqnl, vkb, nkb, & + ntype, natom, rat, adot, & + mxdtyp, mxdatm, mxdgve, mxdlqp, mxddim, mxdbnd) + + do i = 1,3 + do j = 1,3 + sunsym(j,i) = sunsym(j,i) + strnlkb(j,i) + enddo + enddo + + do i = 1,ntype + do j = 1,natom(i) + do k = 1,3 + funsym(k,j,i) = funsym(k,j,i) + fnlkb(k,j,i) + enddo + enddo + enddo + + enddo + + call for_str_sym(funsym, sunsym, fsym, ssym, & + ntype, natom, rat, & + ntrans, mtrx, tnp, & + mxdtyp, mxdatm) + +! sums up contributions + + do i = 1,3 + do j = 1,3 + stress(i,j) = strew(i,j) + strxc(i,j) + do k=1,3 + do l=1,3 + stress(i,j) = stress(i,j) + adotm1(i,k) * & + (strhl(k,l) + ssym(k,l)) * adotm1(l,j) + enddo + enddo + enddo + enddo + + do i = 1,ntype + do j = 1,natom(i) + do k = 1,3 + force(k,j,i) = forcew(k,j,i) + do l = 1,3 + force(k,j,i) = force(k,j,i) + adotm1(k,l) & + * (fsym(l,j,i) + floc(l,j,i)) + enddo + enddo + enddo + enddo + + deallocate(occp) + + deallocate(floc) + deallocate(funsym) + deallocate(fsym) + deallocate(fnlkb) + + return + +end subroutine force_stress_kb diff --git a/Src/cpw_post_process.f90 b/Src/cpw_post_process.f90 index 7421c6a..4d4c177 100644 --- a/Src/cpw_post_process.f90 +++ b/Src/cpw_post_process.f90 @@ -41,7 +41,7 @@ program cpw_post_process ! Driver program version - vdriv = '5.10' + vdriv = '5.11' call tpage(vdriv) diff --git a/Src/read_print/print_energy.f90 b/Src/read_print/print_energy.f90 index dd5da99..24ee253 100644 --- a/Src/read_print/print_energy.f90 +++ b/Src/read_print/print_energy.f90 @@ -11,133 +11,139 @@ ! https://github.com/jlm785/cpw2000 ! !------------------------------------------------------------! -!> prints several energy terms +!> prints several energy terms, forces and stresses +!> +!> \author Jose Luis Martins +!> \version 5.10 +!> \date 15 january 1999, 21 February 2024. +!> \copyright GNU Public License v2 - subroutine print_energy(ipr,entype,energy,force,stress, & - & adot,ntype,natom,nameat, & - & mxdtyp,mxdatm) + subroutine print_energy(ipr, entype, energy, force, stress, & + adot, ntype, natom, nameat, & + mxdtyp, mxdatm) -! written 15 january 1999. jlm -! modified for f90, 21 October 2015. JLM -! Modified, documentation, June 2019. JLM -! copyright INESC-MN/Jose Luis Martins +! written 15 January 1999. jlm +! modified for f90, 21 October 2015. JLM +! Modified, documentation, June 2019. JLM +! Modified, ipr, indentation, 21 February 2024. JLM -! version 4.94 of cpw -! version 1.5 of md + implicit none - implicit none - - integer, parameter :: REAL64 = selected_real_kind(12) - -! input - - integer, intent(in) :: mxdtyp !< array dimension of types of atoms - integer, intent(in) :: mxdatm !< array dimension of number of atoms of a given type - - integer, intent(in) :: ipr !< should be equal to one if information is to be printed. - character(len=*),intent(in) :: entype !< type of energy contribution - real(REAL64), intent(in) :: adot(3,3) !< metric in direct space - integer, intent(in) :: ntype !< number of types of atoms - integer, intent(in) :: natom(mxdtyp) !< number of atoms of type i - character(len=2), intent(in) :: nameat(mxdtyp) !< chemical symbol for the type i - - real(REAL64), intent(in) :: energy !< energy - real(REAL64), intent(in) :: force(3,mxdatm,mxdtyp) !< k-th component (in contravariant lattice coordinates) of the force of the n-th atom of type i - real(REAL64), intent(in) :: stress(3,3) !< stress tensor (in contravariant lattice coordinates) - - -! local variables - - real(REAL64) :: strcar(3,3) - real(REAL64) :: avec(3,3),bvec(3,3),bdot(3,3),vcell - real(REAL64) :: press - real(REAL64) :: car(3) - -! parameters - - real(REAL64), parameter :: ZERO = 0.0_REAL64 - real(REAL64), parameter :: AUTOGPA = 29421.58_REAL64 - -! counters - - integer :: i,j,i1,i2,j1,j2 - integer :: nt,ntt,ja,jmax - - - if(ipr /= 1) return - -! writes energy - - write(6,*) - write(6,*) ' Energy, force, stress from ',entype - write(6,*) - write(6,'(5x,f15.8,5x," energy ",a50)') energy,entype - -! writes stress tensor - - call adot_to_avec_sym(adot,avec,bvec) - call adot_to_bdot(adot,vcell,bdot) - - do i1=1,3 - do i2=1,3 - strcar(i1,i2) = ZERO - do j1=1,3 - do j2=1,3 - strcar(i1,i2) = strcar(i1,i2) + & - & avec(i1,j1) * stress(j1,j2) * avec(i2,j2) - enddo - enddo - strcar(i1,i2) = strcar(i1,i2) * autogpa / vcell - enddo - enddo - - write(6,*) - write(6,'(8x,"Contravariant stress tensor (a.u.)",15x, & - & "Cartesian stress (GPa)")') - do i=1,3 - write (6,'(5x,3(2x,f10.6),5x,3(2x,e14.6)," stress ",i1,a20)') & - & (stress(i,j),j=1,3),(strcar(i,j),j=1,3),i,entype - enddo - write(6,*) - -! pressure - - press = zero - do i=1,3 - do j=1,3 - press = press + adot(i,j)*stress(j,i) - enddo - enddo - press = press/(3*vcell) - write(6,'(5x,f15.8,5x,f15.8,5x," pressure (au and GPa) ",a20)') & - & press,press*autogpa,entype - -! force - - write(6,*) - write(6,'(10x,"Force (Lattice coord.)",12x, & - & "Force (Cartesian coord. a.u)",7x,"no. type ")') - write(6,*) - ntt = 0 - do nt=1,ntype - jmax = natom(nt) - do ja=1,jmax - ntt = ntt + 1 - car(1) = avec(1,1)*force(1,ja,nt) + & - & avec(1,2)*force(2,ja,nt) + & - & avec(1,3)*force(3,ja,nt) - car(2) = avec(2,1)*force(1,ja,nt) + & - & avec(2,2)*force(2,ja,nt) + & - & avec(2,3)*force(3,ja,nt) - car(3) = avec(3,1)*force(1,ja,nt) + & - & avec(3,2)*force(2,ja,nt) + & - & avec(3,3)*force(3,ja,nt) - write(6,'(2x,3(2x,f9.5),3x,3(1x,e12.5),1x,i3,3x,a2,3x, & - & " force ",a20)') (force(i,ja,nt),i=1,3), & - & (car(i),i=1,3),ntt,nameat(nt),entype - enddo - enddo - - return - end subroutine print_energy + integer, parameter :: REAL64 = selected_real_kind(12) + +! input + + integer, intent(in) :: mxdtyp !< array dimension of types of atoms + integer, intent(in) :: mxdatm !< array dimension of number of atoms of a given type + + integer, intent(in) :: ipr !< should be equal to one if information is to be printed. + character(len=*),intent(in) :: entype !< type of energy contribution + real(REAL64), intent(in) :: adot(3,3) !< metric in direct space + integer, intent(in) :: ntype !< number of types of atoms + integer, intent(in) :: natom(mxdtyp) !< number of atoms of type i + character(len=2), intent(in) :: nameat(mxdtyp) !< chemical symbol for the type i + + real(REAL64), intent(in) :: energy !< energy + real(REAL64), intent(in) :: force(3,mxdatm,mxdtyp) !< k-th component (in contravariant lattice coordinates) of the force of the n-th atom of type i + real(REAL64), intent(in) :: stress(3,3) !< stress tensor (in contravariant lattice coordinates) + + +! local variables + + real(REAL64) :: strcar(3,3) + real(REAL64) :: avec(3,3),bvec(3,3),bdot(3,3),vcell + real(REAL64) :: press + real(REAL64) :: car(3) + +! parameters + + real(REAL64), parameter :: ZERO = 0.0_REAL64 + real(REAL64), parameter :: AUTOGPA = 29421.58_REAL64 + +! counters + + integer :: i, j, i1, i2, j1, j2 + integer :: nt, ntt, ja, jmax + + +! keep for compatibility. Eliminate in the future... + + if(ipr < 1) return + +! writes energy + + write(6,*) + write(6,*) ' Energy, force, stress from ', entype + write(6,*) + write(6,'(5x,f15.8,5x," energy ",a50)') energy, entype + +! writes stress tensor + + call adot_to_avec_sym(adot,avec,bvec) + call adot_to_bdot(adot,vcell,bdot) + + do i1 = 1,3 + do i2 = 1,3 + strcar(i1,i2) = ZERO + do j1 = 1,3 + do j2 = 1,3 + strcar(i1,i2) = strcar(i1,i2) + & + avec(i1,j1) * stress(j1,j2) * avec(i2,j2) + enddo + enddo + strcar(i1,i2) = strcar(i1,i2) * autogpa / vcell + enddo + enddo + + write(6,*) + write(6,'(8x,"Contravariant stress tensor (a.u.)",15x, & + & "Cartesian stress (GPa)")') + do i = 1,3 + write (6,'(5x,3(2x,f10.6),5x,3(2x,e14.6)," stress ",i1,a50)') & + (stress(i,j),j=1,3), (strcar(i,j),j=1,3), i, entype + enddo + write(6,*) + +! pressure + + press = zero + do i = 1,3 + do j = 1,3 + press = press + adot(i,j)*stress(j,i) + enddo + enddo + press = press/(3*vcell) + write(6,'(5x,f15.8,5x,f15.8,5x," pressure (au and GPa) ",a20)') & + press, press*autogpa, entype + +! force + + write(6,*) + write(6,'(10x,"Force (Lattice coord.)",12x, & + & "Force (Cartesian coord. a.u)",7x,"no. type ")') + write(6,*) + ntt = 0 + + do nt = 1 ,ntype + jmax = natom(nt) + do ja = 1,jmax + ntt = ntt + 1 + car(1) = avec(1,1)*force(1,ja,nt) + & + avec(1,2)*force(2,ja,nt) + & + avec(1,3)*force(3,ja,nt) + car(2) = avec(2,1)*force(1,ja,nt) + & + avec(2,2)*force(2,ja,nt) + & + avec(2,3)*force(3,ja,nt) + car(3) = avec(3,1)*force(1,ja,nt) + & + avec(3,2)*force(2,ja,nt) + & + avec(3,3)*force(3,ja,nt) + write(6,'(2x,3(2x,f9.5),3x,3(1x,e12.5),1x,i3,3x,a2,3x, & + & " force ",a50)') (force(i,ja,nt),i=1,3), & + (car(i),i=1,3), ntt, nameat(nt), entype + enddo + enddo + + return + +end subroutine print_energy diff --git a/Src/read_print/version.f90 b/Src/read_print/version.f90 index 62c8521..e25e401 100644 --- a/Src/read_print/version.f90 +++ b/Src/read_print/version.f90 @@ -31,7 +31,7 @@ subroutine version(cpwversion, ldevel) character(len=4), intent(out) :: cpwversion !< hardcoded library version logical, intent(out) :: ldevel !< development branch, minor version may be incompatible - cpwversion = '5.10' + cpwversion = '5.11' ldevel = .TRUE. return diff --git a/Tools/pre_relax_vff.f90 b/Tools/pre_relax_vff.f90 index 64e1590..0a7ba27 100644 --- a/Tools/pre_relax_vff.f90 +++ b/Tools/pre_relax_vff.f90 @@ -730,10 +730,10 @@ subroutine vff_constants(iprint, iowrite, & ! Added guessed constants for Pb, May 22, 2020. JLM ! and Arithmetic mean instead of geometric for carbides. JLM ! Added constants from doi:10.1016/j.physe.2009.11.035 +! Bug in mean for carbides (other than SiC), consistency with rede. 21 February 2024. JLM ! copyright INESC-MN/Jose Luis Martins/C.L. Reis -! version 4.99 of pw -! version 1.6.0 of md +! version 5.10 of pw implicit none @@ -767,6 +767,7 @@ subroutine vff_constants(iprint, iowrite, & real(REAL64), parameter :: ANG = 1E-10 real(REAL64), parameter :: UM = 1.0_REAL64 + real(REAL64), parameter :: S3O4 = sqrt(3.0_REAL64) / 4.0_REAL64 ! same species alpha,dist,beta (group IV) @@ -775,15 +776,15 @@ subroutine vff_constants(iprint, iowrite, & nn2 = (n*(n+1))/2 if(nameat(n) == 'C ' .or. nameat(n) == ' C') then alfa(nn2) = 129.33_REAL64 - dist(nn2) = 1.545_REAL64 + dist(nn2) = S3O4*3.5668_REAL64 beta(n,nn2) = 0.655_REAL64*alfa(nn2) elseif(nameat(n) == 'Si') then alfa(nn2) = 48.50_REAL64 - dist(nn2) = 2.352_REAL64 + dist(nn2) = S3O4*5.4310_REAL64 beta(n,nn2) = 0.285_REAL64*alfa(nn2) elseif(nameat(n) == 'Ge') then alfa(nn2) = 38.67_REAL64 - dist(nn2) = 2.450_REAL64 + dist(nn2) = S3O4*5.6579_REAL64 beta(n,nn2) = 0.294_REAL64*alfa(nn2) elseif(nameat(n) == 'Sn') then alfa(nn2) = 25.45_REAL64 @@ -816,25 +817,25 @@ subroutine vff_constants(iprint, iowrite, & & nameat(m) == 'Si') .or. (nameat(n) == 'Si' .and. & & (nameat(m) == 'C ' .or. nameat(m) == ' C'))) then alfa(nm2) = 88.0_REAL64 - dist(nm2) = 1.888_REAL64 + dist(nm2) = S3O4*4.3596_REAL64 beta(n,mm2) = 0.54_REAL64*alfa(nm2) elseif(((nameat(n) == 'C ' .or. nameat(n) == ' C') .and. & & nameat(m) == 'Ge') .or. (nameat(n) == 'Ge' .and. & & (nameat(m) == 'C ' .or. nameat(m) == ' C'))) then alfa(nm2) = (alfa(nn2)+alfa(mm2))/2 - dist(nm2) = (dist(nn2)*dist(mm2))/2 + dist(nm2) = sqrt(dist(nn2)*dist(mm2)) beta(n,mm2) = (beta(n,nn2)+beta(m,mm2))/2 elseif(((nameat(n) == 'C ' .or. nameat(n) == ' C') .and. & & nameat(m) == 'Sn') .or. (nameat(n) == 'Sn' .and. & & (nameat(m) == 'C ' .or. nameat(m) == ' C'))) then alfa(nm2) = (alfa(nn2)+alfa(mm2))/2 - dist(nm2) = (dist(nn2)*dist(mm2))/2 + dist(nm2) = sqrt(dist(nn2)*dist(mm2)) beta(n,mm2) = (beta(n,nn2)+beta(m,mm2))/2 elseif(((nameat(n) == 'C ' .or. nameat(n) == ' C') .and. & & nameat(m) == 'Pb') .or. (nameat(n) == 'Pb' .and. & & (nameat(m) == 'C ' .or. nameat(m) == ' C'))) then alfa(nm2) = (alfa(nn2)+alfa(mm2))/2 - dist(nm2) = (dist(nn2)*dist(mm2))/2 + dist(nm2) = sqrt(dist(nn2)*dist(mm2)) beta(n,mm2) = (beta(n,nn2)+beta(m,mm2))/2 elseif((nameat(n) == 'Si' .and. nameat(m) == 'Ge') .or. & & (nameat(n) == 'Ge' .and. nameat(m) == 'Si')) then @@ -907,7 +908,7 @@ subroutine vff_constants(iprint, iowrite, & elseif((nameat(n) == 'Ga' .and. nameat(m) == 'As') .or. & & (nameat(n) == 'As' .and. nameat(m) == 'Ga')) then alfa(nm2) = 41.19_REAL64 - dist(nm2) = 2.448_REAL64 + dist(nm2) = S3O4*5.6532_REAL64 beta(n,mm2) = 0.217_REAL64*alfa(nm2) elseif((nameat(n) == 'Ga' .and. nameat(m) == 'Sb') .or. & & (nameat(n) == 'Sb' .and. nameat(m) == 'Ga')) then @@ -925,7 +926,7 @@ subroutine vff_constants(iprint, iowrite, & & nameat(m) == ' P')) .or. ((nameat(n) == 'P ' .or. & & nameat(n) == ' P') .and. nameat(m) == 'In')) then alfa(nm2) = 43.04_REAL64 - dist(nm2) = 2.541_REAL64 + dist(nm2) = S3O4*5.8687_REAL64 beta(n,mm2) = 0.145_REAL64*alfa(nm2) elseif((nameat(n) == 'In' .and. nameat(m) == 'As') .or. & & (nameat(n) == 'As' .and. nameat(m) == 'In')) then