diff --git a/src/dftbp/dftb/repulsive/chimesrep.F90 b/src/dftbp/dftb/repulsive/chimesrep.F90 index 1b7b45f9df..ca5e2675f4 100644 --- a/src/dftbp/dftb/repulsive/chimesrep.F90 +++ b/src/dftbp/dftb/repulsive/chimesrep.F90 @@ -7,7 +7,7 @@ #:include 'common.fypp' -!> Implements a repulsive correction using the ChIMES force field +!> Implements a repulsive correction using the ChIMES force field. module dftbp_dftb_repulsive_chimesrep use dftbp_common_accuracy, only : dp use dftbp_common_constants, only : AA__Bohr, Bohr__AA, kcal_mol__Hartree, Hartree__kcal_mol @@ -22,20 +22,27 @@ module dftbp_dftb_repulsive_chimesrep #:if WITH_CHIMES - !> Input for the Chimes repulsive calcualtor + !> Input for the ChIMES repulsive calculator. type :: TChimesRepInp - !> Name of the Chimes input file - character(:), allocatable :: chimesFile + !> Name of the ChIMES input file + character(len=:), allocatable :: chimesFile + + !> Real-space box to enable ChIMES calculations for non-periodic systems. In this case, a + !! user-specified padding is parsed, in order to evaluate the molecule in a large box. + !! The user is responsible for ensuring that this padding is greater than the largest occurring + !! ChIMES cutoff, so that nothing leaks out (in the future, it should be possible to at least + !! query the maximum cutoff via the ChIMES interface, so that this step can be automated). + real(dp), allocatable :: box(:) end type TChimesRepInp - !> Contains the necessary data to query the ChiIMES force field + !> Contains the necessary data to query the ChIMES force field type, extends(TRepulsive) :: TChimesRep private - !> Chimes calculator + !> ChIMES calculator type(TChimesCalc) :: chimesCalc !> Lattice vector buffer @@ -77,22 +84,22 @@ module dftbp_dftb_repulsive_chimesrep #:if WITH_CHIMES - !> Initialises ChIMES repulsive + !> Initialises ChIMES repulsive. subroutine TChimesRep_init(this, input, speciesNames, species0) !> Instance on exit type(TChimesRep), intent(out) :: this - !> Name of the Chimes input file + !> Name of the ChIMES input file type(TChimesRepInp), intent(in) :: input !> Name of the species in the system. Shape: [nSpecies] - character(*), intent(in) :: speciesNames(:) + character(len=*), intent(in) :: speciesNames(:) !> Species index of each atom in the central cell. Shape: [nAtom] integer, intent(in) :: species0(:) - integer :: nAtom + integer :: nAtom, ii call TChimesCalc_init(this%chimesCalc, input%chimesFile, 1, 0) nAtom = size(species0) @@ -100,13 +107,22 @@ subroutine TChimesRep_init(this, input, speciesNames, species0) allocate(this%gradients(3, nAtom), source=0.0_dp) call this%updateSpecies(speciesNames, species0) + ! Non-periodic workaround + if (allocated(input%box)) then + ! Lattice vectors/padding in Bohr, converted to Angström when handed over to ChIMES + this%latVecs(:,:) = 0.0_dp + do ii = 1, 3 + this%latVecs(ii, ii) = input%box(ii) + end do + end if + end subroutine TChimesRep_init !> Notifies the calculator, that order of the species has been changed. subroutine TChimesRep_updateSpecies(this, speciesNames, species) - !> Instance. + !> Instance class(TChimesRep), intent(inout) :: this !> Name of the species in the system. Shape: [nSpecies] @@ -120,7 +136,7 @@ subroutine TChimesRep_updateSpecies(this, speciesNames, species) end subroutine TChimesRep_updateSpecies - !> Returns the real space cutoff (actually 0, as Chimes does not make use of the neighbour list) + !> Returns the real space cutoff (actually 0, as ChIMES does not make use of the neighbour list). function TChimesRep_getRCutOff(this) result(cutOff) !> Instance @@ -163,7 +179,7 @@ subroutine TChimesRep_updateCoords(this, coords, species, img2CentCell, neighbou !> Mapping of atoms into the central cell. Shape: [nAllAtom] integer, intent(in) :: img2CentCell(:) - !> Neighbour list. + !> Neighbour list type(TNeighbourList), intent(in) :: neighbourList real(dp) :: energy @@ -176,9 +192,9 @@ subroutine TChimesRep_updateCoords(this, coords, species, img2CentCell, neighbou call this%chimesCalc%calculate(coords(:, 1:nAtom) * Bohr__AA, this%latVecs * Bohr__AA,& & energy, this%gradients, this%stress) energy = energy * kcal_mol__Hartree - ! Chimes only returns total energy, distribute it equally on each atom. + ! ChIMES only returns total energy, distribute it equally on each atom. this%energyPerAtom(:) = energy / real(size(this%energyPerAtom), dp) - ! Chimes returns forces, not gradients + ! ChIMES returns forces, not gradients this%gradients(:,:) = -this%gradients(:,:) * (kcal_mol__Hartree / AA__Bohr) this%stress(:,:) = this%stress * (kcal_mol__Hartree / AA__Bohr**3) @@ -189,7 +205,7 @@ end subroutine TChimesRep_updateCoords subroutine TChimesRep_getEnergy(this, coords, species, img2CentCell, neighbourList, Eatom,& & Etotal, iAtInCentralRegion) - !> Instance. + !> Instance class(TChimesRep), intent(in) :: this !> All atomic coordinates @@ -224,19 +240,19 @@ subroutine TChimesRep_getEnergy(this, coords, species, img2CentCell, neighbourLi end subroutine TChimesRep_getEnergy - !> Returns the gradients of the repulsive interaction + !> Returns the gradients of the repulsive interaction. subroutine TChimesRep_getGradients(this, coords, species, img2CentCell, neighbourList, grads) !> Instance class(TChimesRep), intent(in) :: this - !> coordinates (x,y,z, all atoms including possible images) + !> Coordinates (x,y,z, all atoms including possible images) real(dp), intent(in) :: coords(:,:) - !> Species of atoms in the central cell. + !> Species of atoms in the central cell integer, intent(in) :: species(:) - !> Index of each atom in the central cell which the atom is mapped on. + !> Index of each atom in the central cell which the atom is mapped on integer, intent(in) :: img2CentCell(:) !> Neighbour list @@ -250,29 +266,29 @@ subroutine TChimesRep_getGradients(this, coords, species, img2CentCell, neighbou end subroutine TChimesRep_getGradients - !> Returns the stress tensor contribution from the repulsive term + !> Returns the stress tensor contribution from the repulsive term. subroutine TChimesRep_getStress(this, coords, species, img2CentCell, neighbourList, cellVol,& & stress) !> Instance class(TChimesRep), intent(in) :: this - !> coordinates (x,y,z, all atoms including possible images) + !> Coordinates (x,y,z, all atoms including possible images) real(dp), intent(in) :: coords(:,:) - !> Species of atoms in the central cell. + !> Species of atoms in the central cell integer, intent(in) :: species(:) - !> indexing array for periodic image atoms + !> Indexing array for periodic image atoms integer, intent(in) :: img2CentCell(:) !> Neighbour list type(TNeighbourList), intent(in) :: neighbourList - !> cell volume. + !> Cell volume real(dp), intent(in) :: cellVol - !> stress tensor + !> Stress tensor real(dp), intent(out) :: stress(:,:) stress(:,:) = this%stress @@ -281,7 +297,7 @@ end subroutine TChimesRep_getStress #:else - !> Dummy initializer in case code was compiled without ChIMES + !> Dummy initializer in case code was compiled without ChIMES. subroutine TChimesRep_init() end subroutine TChimesRep_init diff --git a/src/dftbp/dftbplus/initprogram.F90 b/src/dftbp/dftbplus/initprogram.F90 index b875599e68..8493e4d0fc 100644 --- a/src/dftbp/dftbplus/initprogram.F90 +++ b/src/dftbp/dftbplus/initprogram.F90 @@ -6360,18 +6360,15 @@ subroutine initRepulsive_(nAtom, isPeriodic, isHelical, pairRepulsives, chimesIn @:CREATE_CLASS(repulsive, TTwoBodyRep, TTwoBodyRep_init, twoBodyInp) call repulsiveList%push(repulsive) - #:if WITH_CHIMES - if (allocated(chimesInp)) then - if (.not. isPeriodic) then - call error("ChIMES repulsives currently require periodic boundary conditions") - end if - if (isHelical) then - call error("ChIMES repulsive is not compatible with helical boundary conditions") - end if - @:CREATE_CLASS(repulsive, TChimesRep, TChimesRep_init, chimesInp, speciesNames, species0) - call repulsiveList%push(repulsive) + #:if WITH_CHIMES + if (allocated(chimesInp)) then + if (isHelical) then + call error("ChIMES repulsive is not compatible with helical boundary conditions.") end if - #:endif + @:CREATE_CLASS(repulsive, TChimesRep, TChimesRep_init, chimesInp, speciesNames, species0) + call repulsiveList%push(repulsive) + end if + #:endif ! If multiple repulsives, wrap via container, otherwise use the one directly if (repulsiveList%size() > 1) then diff --git a/src/dftbp/dftbplus/parser.F90 b/src/dftbp/dftbplus/parser.F90 index d95c93c4b1..4b1a0f0b94 100644 --- a/src/dftbp/dftbplus/parser.F90 +++ b/src/dftbp/dftbplus/parser.F90 @@ -1434,7 +1434,7 @@ subroutine readDFTBHam(node, ctrl, geo, slako, poisson) end if end select - call parseChimes(node, ctrl%chimesRepInput) + call parseChimes(node, geo, ctrl%chimesRepInput) ! SCC call getChildValue(node, "SCC", ctrl%tSCC, .false.) @@ -8088,35 +8088,50 @@ subroutine readSpinTuning(node, ctrl, nType) end subroutine readSpinTuning - !> Parses Chimes related options. - subroutine parseChimes(root, chimesRepInput) - type(fnode), pointer, intent(in) :: root - type(TChimesRepInp), allocatable, intent(out) :: chimesRepInput + !> Parses ChIMES related options. + subroutine parseChimes(node, geo, chimesRepInput) + + !> Node to get the information from + type(fnode), intent(in), pointer :: node + + !> Geometry structure + type(TGeometry), intent(in) :: geo + + !> Input structure for the ChIMES repulsive calculator + type(TChimesRepInp), intent(out), allocatable :: chimesRepInput type(fnode), pointer :: chimes - type(string) :: buffer #:if WITH_CHIMES + type(string) :: buffer, modifier type(string), allocatable :: searchPath(:) - #:endif - character(len=:), allocatable :: chimesFile + type(fnode), pointer :: child + real(dp) :: padding + #:endif - call getChild(root, "Chimes", chimes, requested=.false.) + call getChild(node, "Chimes", chimes, requested=.false.) if (.not. associated(chimes)) return - #:if WITH_CHIMES - allocate(chimesRepInput) - call getChildValue(chimes, "ParameterFile", buffer, default="chimes.dat") - chimesFile = unquote(char(buffer)) - call getParamSearchPath(searchPath) - call findFile(searchPath, chimesFile, chimesRepInput%chimesFile) - if (.not. allocated(chimesRepInput%chimesFile)) then - call error("Could not find ChIMES parameter file '" // chimesFile // "'") - end if - #:else - call detailedError(chimes, "ChIMES repuslive correction requested, but code was compiled& - & without ChIMES support") - #:endif + + #:if WITH_CHIMES + allocate(chimesRepInput) + call getChildValue(chimes, "ParameterFile", buffer, default="chimes.dat") + chimesFile = unquote(char(buffer)) + call getParamSearchPath(searchPath) + call findFile(searchPath, chimesFile, chimesRepInput%chimesFile) + if (.not. allocated(chimesRepInput%chimesFile)) then + call error("Could not find ChIMES parameter file '" // chimesFile // "'.") + end if + if (.not. geo%tPeriodic) then + call getChildValue(chimes, "Padding", padding, modifier=modifier, child=child) + call convertUnitHsd(char(modifier), lengthUnits, child, padding) + allocate(chimesRepInput%box(3)) + chimesRepInput%box(:) = maxval(geo%coords, dim=2) - minval(geo%coords, dim=2) + padding + end if + #:else + call detailedError(chimes, "ChIMES repulsive correction requested, but code was compiled& + & without ChIMES support.") + #:endif end subroutine parseChimes