diff --git a/include/xtb.h b/include/xtb.h index 383c264e1..4c11b8cc8 100644 --- a/include/xtb.h +++ b/include/xtb.h @@ -19,10 +19,10 @@ #define XTB_API_ENTRY #define XTB_API_CALL -#define XTB_API_SUFFIX__VERSION_1_0_0 +#define XTB_API_SUFFIX__VERSION_2_0_0 /// Define proprocessor to allow to check for specific API features -#define XTB_API_VERSION 10000 +#define XTB_API_VERSION 20000 #define XTB_VERSION_6_3_0 1 #define XTB_VERSION_6_3_1 1 #define XTB_VERSION_6_3_2 1 @@ -69,7 +69,7 @@ typedef struct _xtb_TResults* xtb_TResults; /// Returns API version as 10000 * major + 100 * minor + 1 * patch extern XTB_API_ENTRY int XTB_API_CALL -xtb_getAPIVersion() XTB_API_SUFFIX__VERSION_1_0_0; +xtb_getAPIVersion() XTB_API_SUFFIX__VERSION_2_0_0; /* * Calculation environment @@ -77,40 +77,40 @@ xtb_getAPIVersion() XTB_API_SUFFIX__VERSION_1_0_0; /// Create new xtb calculation environment object extern XTB_API_ENTRY xtb_TEnvironment XTB_API_CALL -xtb_newEnvironment(void) XTB_API_SUFFIX__VERSION_1_0_0; +xtb_newEnvironment(void) XTB_API_SUFFIX__VERSION_2_0_0; /// Delete a xtb calculation environment object extern XTB_API_ENTRY void XTB_API_CALL -xtb_delEnvironment(xtb_TEnvironment* /* env */) XTB_API_SUFFIX__VERSION_1_0_0; +xtb_delEnvironment(xtb_TEnvironment* /* env */) XTB_API_SUFFIX__VERSION_2_0_0; /// Check current status of calculation environment extern XTB_API_ENTRY int XTB_API_CALL -xtb_checkEnvironment(xtb_TEnvironment /* env */) XTB_API_SUFFIX__VERSION_1_0_0; +xtb_checkEnvironment(xtb_TEnvironment /* env */) XTB_API_SUFFIX__VERSION_2_0_0; /// Show and empty error stack extern XTB_API_ENTRY void XTB_API_CALL xtb_showEnvironment(xtb_TEnvironment /* env */, - const char* /* message */) XTB_API_SUFFIX__VERSION_1_0_0; + const char* /* message */) XTB_API_SUFFIX__VERSION_2_0_0; /// Return and empty error stack extern XTB_API_ENTRY void XTB_API_CALL xtb_getError(xtb_TEnvironment /* env */, char* /* buffer */, - const int* /* buffersize */) XTB_API_SUFFIX__VERSION_1_0_0; + const int* /* buffersize */) XTB_API_SUFFIX__VERSION_2_0_0; /// Bind output from this environment extern XTB_API_ENTRY void XTB_API_CALL xtb_setOutput(xtb_TEnvironment /* env */, - const char* /* filename */) XTB_API_SUFFIX__VERSION_1_0_0; + const char* /* filename */) XTB_API_SUFFIX__VERSION_2_0_0; /// Release output unit from this environment extern XTB_API_ENTRY void XTB_API_CALL -xtb_releaseOutput(xtb_TEnvironment /* env */) XTB_API_SUFFIX__VERSION_1_0_0; +xtb_releaseOutput(xtb_TEnvironment /* env */) XTB_API_SUFFIX__VERSION_2_0_0; /// Set verbosity of calculation output extern XTB_API_ENTRY void XTB_API_CALL xtb_setVerbosity(xtb_TEnvironment /* env */, - int /* verbosity */) XTB_API_SUFFIX__VERSION_1_0_0; + int /* verbosity */) XTB_API_SUFFIX__VERSION_2_0_0; /* * Molecular structure data class @@ -125,18 +125,18 @@ xtb_newMolecule(xtb_TEnvironment /* env */, const double* /* charge in e */, const int* /* uhf */, const double* /* lattice [3][3] */, - const bool* /* periodic [3] */) XTB_API_SUFFIX__VERSION_1_0_0; + const bool* /* periodic [3] */) XTB_API_SUFFIX__VERSION_2_0_0; /// Delete molecular structure data extern XTB_API_ENTRY void XTB_API_CALL -xtb_delMolecule(xtb_TMolecule* /* mol */) XTB_API_SUFFIX__VERSION_1_0_0; +xtb_delMolecule(xtb_TMolecule* /* mol */) XTB_API_SUFFIX__VERSION_2_0_0; /// Update coordinates and lattice parameters (quantities in Bohr) extern XTB_API_ENTRY void XTB_API_CALL xtb_updateMolecule(xtb_TEnvironment /* env */, xtb_TMolecule /* mol */, const double* /* positions [natoms][3] */, - const double* /* lattice [3][3] */) XTB_API_SUFFIX__VERSION_1_0_0; + const double* /* lattice [3][3] */) XTB_API_SUFFIX__VERSION_2_0_0; /* * Singlepoint calculator @@ -144,39 +144,39 @@ xtb_updateMolecule(xtb_TEnvironment /* env */, /// Create new calculator object extern XTB_API_ENTRY xtb_TCalculator XTB_API_CALL -xtb_newCalculator(void) XTB_API_SUFFIX__VERSION_1_0_0; +xtb_newCalculator(void) XTB_API_SUFFIX__VERSION_2_0_0; /// Delete calculator object extern XTB_API_ENTRY void XTB_API_CALL -xtb_delCalculator(xtb_TCalculator* /* calc */) XTB_API_SUFFIX__VERSION_1_0_0; +xtb_delCalculator(xtb_TCalculator* /* calc */) XTB_API_SUFFIX__VERSION_2_0_0; /// Load GFN0-xTB calculator extern XTB_API_ENTRY void XTB_API_CALL xtb_loadGFN0xTB(xtb_TEnvironment /* env */, xtb_TMolecule /* mol */, xtb_TCalculator /* calc */, - char* /* filename */) XTB_API_SUFFIX__VERSION_1_0_0; + char* /* filename */) XTB_API_SUFFIX__VERSION_2_0_0; /// Load GFN1-xTB calculator extern XTB_API_ENTRY void XTB_API_CALL xtb_loadGFN1xTB(xtb_TEnvironment /* env */, xtb_TMolecule /* mol */, xtb_TCalculator /* calc */, - char* /* filename */) XTB_API_SUFFIX__VERSION_1_0_0; + char* /* filename */) XTB_API_SUFFIX__VERSION_2_0_0; /// Load GFN2-xTB calculator extern XTB_API_ENTRY void XTB_API_CALL xtb_loadGFN2xTB(xtb_TEnvironment /* env */, xtb_TMolecule /* mol */, xtb_TCalculator /* calc */, - char* /* filename */) XTB_API_SUFFIX__VERSION_1_0_0; + char* /* filename */) XTB_API_SUFFIX__VERSION_2_0_0; /// Load GFN-FF calculator extern XTB_API_ENTRY void XTB_API_CALL xtb_loadGFNFF(xtb_TEnvironment /* env */, xtb_TMolecule /* mol */, xtb_TCalculator /* calc */, - char* /* filename */) XTB_API_SUFFIX__VERSION_1_0_0; + char* /* filename */) XTB_API_SUFFIX__VERSION_2_0_0; /// Add a solvation model to calculator (requires loaded parametrisation) extern XTB_API_ENTRY void XTB_API_CALL @@ -185,12 +185,13 @@ xtb_setSolvent(xtb_TEnvironment /* env */, char* /* solvent */, int* /* state */, double* /* temp */, - int* /* grid */) XTB_API_SUFFIX__VERSION_1_0_0; + int* /* grid */, + char* /* solvent model */) XTB_API_SUFFIX__VERSION_2_0_0; /// Unset the solvation model extern XTB_API_ENTRY void XTB_API_CALL xtb_releaseSolvent(xtb_TEnvironment /* env */, - xtb_TCalculator /* calc */) XTB_API_SUFFIX__VERSION_1_0_0; + xtb_TCalculator /* calc */) XTB_API_SUFFIX__VERSION_2_0_0; /// Add a external charge potential to calculator (only supported in GFN1/2-xTB) extern XTB_API_ENTRY void XTB_API_CALL @@ -199,37 +200,44 @@ xtb_setExternalCharges(xtb_TEnvironment /* env */, int* /* n */, int* /* numbers [n] */, double* /* charges [n] */, - double* /* positions [n][3] */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* positions [n][3] */) XTB_API_SUFFIX__VERSION_2_0_0; /// Unset the external charge potential extern XTB_API_ENTRY void XTB_API_CALL xtb_releaseExternalCharges(xtb_TEnvironment /* env */, - xtb_TCalculator /* calc */) XTB_API_SUFFIX__VERSION_1_0_0; + xtb_TCalculator /* calc */) XTB_API_SUFFIX__VERSION_2_0_0; /// Set numerical accuracy of calculator in the range of 1000 to 0.0001 extern XTB_API_ENTRY void XTB_API_CALL xtb_setAccuracy(xtb_TEnvironment /* env */, xtb_TCalculator /* calc */, - double /* accuracy */) XTB_API_SUFFIX__VERSION_1_0_0; + double /* accuracy */) XTB_API_SUFFIX__VERSION_2_0_0; /// Set maximum number of iterations for self-consistent TB calculators extern XTB_API_ENTRY void XTB_API_CALL xtb_setMaxIter(xtb_TEnvironment /* env */, xtb_TCalculator /* calc */, - int /* iterations */) XTB_API_SUFFIX__VERSION_1_0_0; + int /* iterations */) XTB_API_SUFFIX__VERSION_2_0_0; /// Set electronic temperature for level filling in tight binding calculators in K extern XTB_API_ENTRY void XTB_API_CALL xtb_setElectronicTemp(xtb_TEnvironment /* env */, xtb_TCalculator /* calc */, - double /* temperature */) XTB_API_SUFFIX__VERSION_1_0_0; + double /* temperature */) XTB_API_SUFFIX__VERSION_2_0_0; + +//// Calculate CPCM-X solvation energy +extern XTB_API_ENTRY void XTB_API_CALL +xtb_cpcmx_calc(xtb_TEnvironment /* env */, + xtb_TMolecule /* mol */, + xtb_TCalculator /* calc */, + xtb_TResults /* res */) XTB_API_SUFFIX__VERSION_2_0_0; /// Perform singlepoint calculation extern XTB_API_ENTRY void XTB_API_CALL xtb_singlepoint(xtb_TEnvironment /* env */, xtb_TMolecule /* mol */, xtb_TCalculator /* calc */, - xtb_TResults /* res */) XTB_API_SUFFIX__VERSION_1_0_0; + xtb_TResults /* res */) XTB_API_SUFFIX__VERSION_2_0_0; /* * Calculation results @@ -237,81 +245,87 @@ xtb_singlepoint(xtb_TEnvironment /* env */, /// Create new singlepoint results object extern XTB_API_ENTRY xtb_TResults XTB_API_CALL -xtb_newResults(void) XTB_API_SUFFIX__VERSION_1_0_0; +xtb_newResults(void) XTB_API_SUFFIX__VERSION_2_0_0; /// Delete singlepoint results object extern XTB_API_ENTRY void XTB_API_CALL -xtb_delResults(xtb_TResults* /* res */) XTB_API_SUFFIX__VERSION_1_0_0; +xtb_delResults(xtb_TResults* /* res */) XTB_API_SUFFIX__VERSION_2_0_0; /// Create copy from a singlepoint results object extern XTB_API_ENTRY xtb_TResults XTB_API_CALL -xtb_copyResults(xtb_TResults /* res */) XTB_API_SUFFIX__VERSION_1_0_0; +xtb_copyResults(xtb_TResults /* res */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for energy in Hartree extern XTB_API_ENTRY void XTB_API_CALL xtb_getEnergy(xtb_TEnvironment /* env */, xtb_TResults /* res */, - double* /* energy */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* energy */) XTB_API_SUFFIX__VERSION_2_0_0; + +/// Query singlepoint results object for solvation energy in Hartree +extern XTB_API_ENTRY void XTB_API_CALL +xtb_getSolvationEnergy(xtb_TEnvironment /* env */, + xtb_TResults /* res */, + double* /* energy */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for gradient in Hartree / Bohr extern XTB_API_ENTRY void XTB_API_CALL xtb_getGradient(xtb_TEnvironment /* env */, xtb_TResults /* res */, - double* /* gradient [natoms][3] */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* gradient [natoms][3] */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for pc gradient in Hartree / Bohr extern XTB_API_ENTRY void XTB_API_CALL xtb_getPCGradient(xtb_TEnvironment /* env */, xtb_TResults /* res */, - double* /* gradient [natoms][3] */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* gradient [natoms][3] */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for virial in Hartree extern XTB_API_ENTRY void XTB_API_CALL xtb_getVirial(xtb_TEnvironment /* env */, xtb_TResults /* res */, - double* /* virial [3][3] */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* virial [3][3] */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for dipole in e Bohr extern XTB_API_ENTRY void XTB_API_CALL xtb_getDipole(xtb_TEnvironment /* env */, xtb_TResults /* res */, - double* /* dipole [3] */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* dipole [3] */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for partial charges in e extern XTB_API_ENTRY void XTB_API_CALL xtb_getCharges(xtb_TEnvironment /* env */, xtb_TResults /* res */, - double* /* charges [natoms] */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* charges [natoms] */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for bond orders extern XTB_API_ENTRY void XTB_API_CALL xtb_getBondOrders(xtb_TEnvironment /* env */, xtb_TResults /* res */, - double* /* wbo [natoms][natoms] */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* wbo [natoms][natoms] */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for the number of basis functions extern XTB_API_ENTRY void XTB_API_CALL xtb_getNao(xtb_TEnvironment /* env */, xtb_TResults /* res */, - int* /* nao */) XTB_API_SUFFIX__VERSION_1_0_0; + int* /* nao */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for orbital energies in Hartree [nao] extern XTB_API_ENTRY void XTB_API_CALL xtb_getOrbitalEigenvalues(xtb_TEnvironment /* env */, xtb_TResults /* res */, - double* /* emo */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* emo */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for occupation numbers [nao] extern XTB_API_ENTRY void XTB_API_CALL xtb_getOrbitalOccupations(xtb_TEnvironment /* env */, xtb_TResults /* res */, - double* /* focc */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* focc */) XTB_API_SUFFIX__VERSION_2_0_0; /// Query singlepoint results object for orbital coefficients [nao][nao] extern XTB_API_ENTRY void XTB_API_CALL xtb_getOrbitalCoefficients(xtb_TEnvironment /* env */, xtb_TResults /* res */, - double* /* c */) XTB_API_SUFFIX__VERSION_1_0_0; + double* /* c */) XTB_API_SUFFIX__VERSION_2_0_0; #ifdef __cplusplus } diff --git a/src/api/calculator.f90 b/src/api/calculator.f90 index 1fe9c706a..bb80e8345 100644 --- a/src/api/calculator.f90 +++ b/src/api/calculator.f90 @@ -338,7 +338,7 @@ end subroutine loadGFN2xTB_api !> Add a solvation model to calculator (requires loaded parametrisation) -subroutine setSolvent_api(venv, vcalc, charptr, state, temperature, grid) & +subroutine setSolvent_api(venv, vcalc, charptr, state, temperature, grid, charptr2) & & bind(C, name="xtb_setSolvent") !DEC$ ATTRIBUTES DLLEXPORT :: setSolvent_api character(len=*), parameter :: source = 'xtb_api_setSolvent' @@ -346,11 +346,13 @@ subroutine setSolvent_api(venv, vcalc, charptr, state, temperature, grid) & type(VEnvironment), pointer :: env type(c_ptr), value :: vcalc type(VCalculator), pointer :: calc - character(kind=c_char), intent(in) :: charptr(*) + character(kind=c_char), intent(in) :: charptr(*) ! Solvent integer(c_int), intent(in), optional :: state real(c_double), intent(in), optional :: temperature integer(c_int), intent(in), optional :: grid + character(kind=c_char), intent(in), optional :: charptr2(*) ! Solvent model character(len=:), allocatable :: solvent + character(len=:), allocatable :: solv_model type(TSolvInput) :: input integer :: gsolvstate, nang real(wp) :: temp @@ -391,6 +393,12 @@ subroutine setSolvent_api(venv, vcalc, charptr, state, temperature, grid) & call c_f_character(charptr, solvent) + if (present(charptr2)) then + call c_f_character(charptr2, solv_model) + else + solv_model = 'gbsa' + end if + ! PGI 20.5 cannot use default constructor with deferred-length characters: ! input = TSolvInput(solvent=solvent, temperature=temp, state=gsolvstate, & ! & nang=nang) @@ -398,8 +406,32 @@ subroutine setSolvent_api(venv, vcalc, charptr, state, temperature, grid) & input%temperature = temp input%state = gsolvstate input%nang = nang + + ! Set solvent model input%alpb = .false. - input%kernel = gbKernel%still + input%cosmo = .false. + input%tmcosmo = .false. + input%kernel = gbKernel%p16 + if (solv_model == 'gbsa') then + input%kernel = gbKernel%still + elseif (solv_model == 'alpb') then + input%alpb = .true. + elseif (solv_model == 'cosmo') then + input%cosmo = .true. + elseif (solv_model == 'tmcosmo') then + input%tmcosmo = .true. + elseif (solv_model == 'cpcmx') then + ! CPCM-X does an initial SCF with COSMO and a special solvent + ! before running a second SCF with the actual solvent. + input%cosmo = .true. + input%solvent = 'infinity' + + input%cpxsolvent = solvent + else + call env%ptr%error("Unknown solvation model", source) + return + end if + call addSolvationModel(env%ptr, calc%ptr, input) call env%ptr%check(exitRun) diff --git a/src/api/interface.f90 b/src/api/interface.f90 index 70b1c1b7c..e626893a0 100644 --- a/src/api/interface.f90 +++ b/src/api/interface.f90 @@ -155,5 +155,156 @@ subroutine singlepoint_api(venv, vmol, vcalc, vres) & end subroutine singlepoint_api +subroutine cpcmx_calc_api(venv, vmol, vcalc, vres) & + & bind(C, name="xtb_cpcmx_calc") + !DEC$ ATTRIBUTES DLLEXPORT :: cpcmx_calc_api + use xtb_solv_cpx, only: TCpcmx + use xtb_type_calculator, only: TCalculator + + character(len=*), parameter :: source = 'xtb_api_cpcmx_calc' + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vmol + type(VMolecule), pointer :: mol + type(c_ptr), value :: vcalc + type(VCalculator), pointer :: calc + type(c_ptr), value :: vres + type(VResults), pointer :: res + type(scc_results) :: spRes + + type(TCpcmx) :: cpx + type(VCalculator) :: cpx_calc + real(c_double) :: energy_gas + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vmol)) then + call env%ptr%error("Molecular structure data is not allocated", source) + return + end if + call c_f_pointer(vmol, mol) + + if (.not.c_associated(vcalc)) then + call env%ptr%error("Singlepoint calculator is not allocated", source) + return + end if + call c_f_pointer(vcalc, calc) + + ! Fail early if CPCM-X solvation model is not set + if (allocated(calc%ptr%solvation)) then + if (.not. allocated(calc%ptr%solvation%cpxsolvent)) then + call env%ptr%error("CPCM-X solvent not set", source) + return + end if + else + call env%ptr%error("No solvation input given", source) + return + end if + + ! Fail early if not using xTB + select type(xtb => calc%ptr) + type is (TGFFCalculator) + call env%ptr%error("CPCM-X is not possible with a force field.", source) + return + end select + + if (.not.allocated(calc%ptr)) then + call env%ptr%error("No calculator loaded for single point", & + & source) + return + end if + + if (.not.c_associated(vres)) then + call env%ptr%error("Calculation results are not allocated", source) + return + end if + call c_f_pointer(vres, res) + + ! check cache, automatically invalidate missmatched data + if (allocated(res%chk)) then + select type(xtb => calc%ptr) + type is(TxTBCalculator) + if (res%chk%wfn%n /= mol%ptr%n .or. res%chk%wfn%n /= xtb%basis%n .or. & + & res%chk%wfn%nao /= xtb%basis%nao .or. & + & res%chk%wfn%nshell /= xtb%basis%nshell) then + deallocate(res%chk) + end if + end select + end if + + if (.not.allocated(res%chk)) then + allocate(res%chk) + ! in case of a new wavefunction cache we have to perform an initial guess + select type(xtb => calc%ptr) + type is(TxTBCalculator) + call newWavefunction(env%ptr, mol%ptr, xtb, res%chk) + end select + end if + + if (.not.allocated(res%energy)) then + allocate(res%energy) + end if + + if (.not.allocated(res%egap)) then + allocate(res%egap) + end if + + if (allocated(res%pcgradient)) then + deallocate(res%pcgradient) + end if + + if (allocated(res%gradient)) then + if (any(shape(res%gradient) /= [3, mol%ptr%n])) then + call env%ptr%warning("Shape missmatch in gradient, reallocating", source) + deallocate(res%gradient) + end if + end if + if (.not.allocated(res%gradient)) then + allocate(res%gradient(3, mol%ptr%n)) + end if + + if (allocated(res%sigma)) then + if (any(shape(res%sigma) /= [3, 3])) then + call env%ptr%warning("Shape missmatch in virial, reallocating", source) + deallocate(res%sigma) + end if + end if + if (.not.allocated(res%sigma)) then + allocate(res%sigma(3, 3)) + end if + + ! Initial COSMO singlepoint calculation + call calc%ptr%singlepoint(env%ptr, mol%ptr, res%chk, env%verbosity, .true., & + & res%energy, res%gradient, res%sigma, res%egap, spRes) + + ! CPCM-X calculation + call cpx%setup(env%ptr, calc%ptr%solvation%cpxsolvent) + cpx_calc = calc + deallocate(cpx_calc%ptr%solvation) + + energy_gas = res%energy + call generic_header(env%ptr%unit, "CPCM-X post-SCF solvation evaluation", 49, 10) + call cpx_calc%ptr%singlepoint(env%ptr, mol%ptr, res%chk, env%verbosity, .false., & + & energy_gas, res%gradient, res%sigma, res%egap, spRes) + + + call cpx%calc_solv(env%ptr, calc%ptr%solvation%cpxsolvent, energy_gas, & + & 0.4_wp, 298.15_wp, 500, 0.0001_wp, spRes%e_total) + call cpx%print(.true.) + + res%energy = spRes%e_total + res%solvation_energy = res%energy - energy_gas + res%dipole = spRes%dipole + + ! Zero out the gradient and sigma (not yet implemented for CPCM-X) + res%gradient = 0.0_wp + res%sigma = 0.0_wp + + endif + +end subroutine cpcmx_calc_api + end module xtb_api_interface diff --git a/src/api/results.f90 b/src/api/results.f90 index 74c99bc50..2c904b918 100644 --- a/src/api/results.f90 +++ b/src/api/results.f90 @@ -42,6 +42,7 @@ module xtb_api_results real(wp), allocatable :: sigma(:, :) real(wp), allocatable :: dipole(:) real(wp), allocatable :: egap + real(wp), allocatable :: solvation_energy end type VResults @@ -135,6 +136,38 @@ subroutine getEnergy_api(venv, vres, dptr) & end subroutine getEnergy_api +!> Query singlepoint results object for solvation energy +subroutine getSolvationEnergy_api(venv, vres, dptr) & + & bind(C, name="xtb_getSolvationEnergy") + !DEC$ ATTRIBUTES DLLEXPORT :: getSolvationEnergy_api + character(len=*), parameter :: source = "xtb_api_getSolvationEnergy" + type(c_ptr), value :: venv + type(VEnvironment), pointer :: env + type(c_ptr), value :: vres + type(VResults), pointer :: res + real(c_double), intent(inout) :: dptr + + if (c_associated(venv)) then + call c_f_pointer(venv, env) + call checkGlobalEnv + + if (.not.c_associated(vres)) then + call env%ptr%error("Results object is not allocated", source) + return + end if + call c_f_pointer(vres, res) + + if (.not.allocated(res%solvation_energy)) then + call env%ptr%error("Solvation energy is not available in results", source) + return + end if + + dptr = res%solvation_energy + + end if + +end subroutine getSolvationEnergy_api + !> Query singlepoint results object for gradient subroutine getGradient_api(venv, vres, dptr) & diff --git a/src/api/version.f90 b/src/api/version.f90 index e8d5b1b4b..851fd70e2 100644 --- a/src/api/version.f90 +++ b/src/api/version.f90 @@ -23,7 +23,7 @@ module xtb_api_version public :: getAPIVersion_api - integer(c_int), parameter :: apiMajor = 1 + integer(c_int), parameter :: apiMajor = 2 integer(c_int), parameter :: apiMinor = 0 integer(c_int), parameter :: apiPatch = 0 diff --git a/src/scf_module.F90 b/src/scf_module.F90 index 7a6ebfb37..a897ecb30 100644 --- a/src/scf_module.F90 +++ b/src/scf_module.F90 @@ -851,17 +851,18 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, & & allocated(solvation), basis, res) endif - if (allocated(solvation)) then - select type(solvation) - type is (TCosmo) - call open_file(ich, "xtb.cosmo", 'w') - call solvation%writeCosmoFile(ich, mol%at, mol%sym, mol%xyz, & - & wfn%q, eel + ep + exb + merge(embd, ed + embd, allocated(scD4))) - call close_file(ich) - end select - end if - endif printing + + ! Need to write xtb.cosmo for CPCM-X, so separate this from the printing block + if (allocated(solvation)) then + select type(solvation) + type is (TCosmo) + call open_file(ich, "xtb.cosmo", 'w') + call solvation%writeCosmoFile(ich, mol%at, mol%sym, mol%xyz, & + & wfn%q, eel + ep + exb + merge(embd, ed + embd, allocated(scD4))) + call close_file(ich) + end select + end if !--------------------------! ! Wiberg-Mayer bond orders ! diff --git a/test/api/c_api_example.c b/test/api/c_api_example.c index a4c04ee49..116faaccd 100644 --- a/test/api/c_api_example.c +++ b/test/api/c_api_example.c @@ -22,12 +22,12 @@ static inline bool check_double(double actual, double expected, double tol, if (fabs(expected - actual) < tol) { return true; } - fprintf(stderr, "FAIL: %s: expected %g, got %g\n", msg, expected, actual); + fprintf(stderr, "FAIL: %s: expected %.10f, got %.10f\n", msg, expected, actual); return false; } int testFirst() { - // first molecule + // first molecule, various solvent models double* q; char* buffer; double* wbo; @@ -51,11 +51,16 @@ int testFirst() { xtb_TCalculator calc = NULL; xtb_TResults res = NULL; double energy; + double solv_energy; double dipole[3]; q = (double*) malloc(natoms * sizeof(double)); wbo = (double*) malloc(natsq * sizeof(double)); buffer = (char*) malloc(buffersize *sizeof(char)); char solvent[] = "h2o"; + char gbsa[] = "gbsa"; + char alpb[] = "alpb"; + char cosmo[] = "cosmo"; + char cpcmx[] = "cpcmx"; if (!check(XTB_API_VERSION, xtb_getAPIVersion(), "API version does not match")) goto error; @@ -104,7 +109,7 @@ int testFirst() { if (!check(wbo[9], 2.89823984265213, 1.0e-8, "Bond order does not match")) goto error; - xtb_setSolvent(env, calc, solvent, NULL, NULL, NULL); + xtb_setSolvent(env, calc, solvent, NULL, NULL, NULL, gbsa); if (xtb_checkEnvironment(env)) goto error; @@ -119,15 +124,119 @@ int testFirst() { if (xtb_checkEnvironment(env)) goto error; - if (!check(energy, -8.38393864716134, 1.0e-9, "Energy does not match")) + if (!check(energy, -8.38393864716134, 1.0e-9, "GBSA Energy does not match")) goto error; - if (!check(q[5], 0.06090868805034, 1.0e-8, "Charge does not match")) + if (!check(q[5], 0.06090868805034, 1.0e-8, "GBSA Charge does not match")) goto error; - if (!check(dipole[2], -0.35455233974705, 1.0e-6, "Dipole does not match")) + if (!check(dipole[2], -0.35455233974705, 1.0e-6, "GBSA Dipole does not match")) goto error; - if (!check(wbo[9], +2.89453979224265, 1.0e-8, "Bond order does not match")) + if (!check(wbo[9], +2.89453979224265, 1.0e-8, "GBSA Bond order does not match")) goto error; + // ALPB + xtb_releaseSolvent(env, calc); + xtb_setSolvent(env, calc, solvent, NULL, NULL, NULL, alpb); + if (xtb_checkEnvironment(env)) + goto error; + + xtb_singlepoint(env, mol, calc, res); + if (xtb_checkEnvironment(env)) + goto error; + + xtb_getEnergy(env, res, &energy); + xtb_getCharges(env, res, q); + xtb_getDipole(env, res, dipole); + xtb_getBondOrders(env, res, wbo); + if (xtb_checkEnvironment(env)) + goto error; + + if (!check(energy, -8.384076843892, 1.0e-9, "ALPB Energy does not match")) + goto error; + if (!check(q[5], 0.0644849340, 1.0e-8, "ALPB Charge does not match")) + goto error; + if (!check(dipole[2], -0.3641866008, 1.0e-6, "ALPB Dipole does not match")) + goto error; + if (!check(wbo[9], 2.8932146955, 1.0e-8, "ALPB Bond order does not match")) + goto error; + + // COSMO + //Sensitive to guess, so reset the results + xtb_delete(res); + res = xtb_newResults(); + if (xtb_checkEnvironment(env)) + goto error; + + xtb_loadGFN2xTB(env, mol, calc, NULL); + xtb_setAccuracy(env, calc, 1.0); + xtb_setElectronicTemp(env, calc, 300.0); + xtb_setMaxIter(env, calc, 30); + if (xtb_checkEnvironment(env)) + goto error; + + xtb_releaseSolvent(env, calc); + xtb_setSolvent(env, calc, solvent, NULL, NULL, NULL, cosmo); + if (xtb_checkEnvironment(env)) + goto error; + + xtb_singlepoint(env, mol, calc, res); + if (xtb_checkEnvironment(env)) + goto error; + + xtb_getEnergy(env, res, &energy); + xtb_getCharges(env, res, q); + xtb_getDipole(env, res, dipole); + xtb_getBondOrders(env, res, wbo); + if (xtb_checkEnvironment(env)) + goto error; + + if (!check(energy, -8.385752224008, 1.0e-9, "COSMO Energy does not match")) + goto error; + if (!check(q[5], 0.0597931738, 1.0e-8, "COSMO Charge does not match")) + goto error; + if (!check(dipole[2], -0.3709356947, 1.0e-6, "COSMO Dipole does not match")) + goto error; + if (!check(wbo[9], 2.8940365143, 1.0e-8, "COSMO Bond order does not match")) + goto error; + + + // CPCMX +#if WITH_CPCMX + //Sensitive to guess, so reset the results + xtb_delete(res); + res = xtb_newResults(); + if (xtb_checkEnvironment(env)) + goto error; + + xtb_releaseSolvent(env, calc); + xtb_setSolvent(env, calc, solvent, NULL, NULL, NULL, cpcmx); + if (xtb_checkEnvironment(env)) + goto error; + + xtb_cpcmx_calc(env, mol, calc, res); + if (xtb_checkEnvironment(env)) + goto error; + + xtb_getEnergy(env, res, &energy); + xtb_getSolvationEnergy(env, res, &solv_energy); + xtb_getCharges(env, res, q); + xtb_getDipole(env, res, dipole); + xtb_getBondOrders(env, res, wbo); + if (xtb_checkEnvironment(env)) + goto error; + + if (!check(energy, -8.383520040765, 1.0e-9, "CPCM-X Energy does not match")) + goto error; + if (!check(solv_energy, -0.0010406558065, 1.0e-8, "CPCM-X Solvation Energy does not match")) + goto error; + if (!check(q[5], 0.0518386448, 1.0e-8, "CPCM-X Charge does not match")) + goto error; + if (!check(dipole[2], -0.2983150104, 1.0e-6, "CPCM-X Dipole does not match")) + goto error; + if (!check(wbo[9], 2.8982388543, 1.0e-8, "CPCM-X Bond order does not match")) + goto error; + +#endif + xtb_delete(res); xtb_delete(calc); xtb_delete(mol); diff --git a/test/api/meson.build b/test/api/meson.build index c13064796..3ef3c6e82 100644 --- a/test/api/meson.build +++ b/test/api/meson.build @@ -19,7 +19,10 @@ test( executable( 'xtb_c_test', sources: files('c_api_example.c'), - dependencies: xtb_dep_static + dependencies: xtb_dep_static, + c_args: [ + '-DWITH_CPCMX=@0@'.format(cpx_dep.found() ? 1 : 0), + ], ), env: xtbenv, )