# Interfacing with a Fortran host¶

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 ! ! PCMSolver, an API for the Polarizable Continuum Model ! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. ! ! This file is part of PCMSolver. ! ! PCMSolver is free software: you can redistribute it and/or modify ! it under the terms of the GNU Lesser General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! PCMSolver is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public License ! along with PCMSolver. If not, see . ! ! For information on the complete list of contributors to the ! PCMSolver API, see: ! program pcm_fortran_host use, intrinsic :: iso_c_binding use, intrinsic :: iso_fortran_env, only:output_unit, error_unit use pcmsolver use utilities use testing implicit none type(c_ptr) :: pcm_context integer(c_int) :: nr_nuclei real(c_double), allocatable :: charges(:) real(c_double), allocatable :: coordinates(:) integer(c_int) :: symmetry_info(4) type(PCMInput) :: host_input logical :: log_open, log_exist character(kind=c_char, len=*), parameter :: mep_lbl = 'NucMEP' character(kind=c_char, len=*), parameter :: asc_lbl = 'NucASC' character(kind=c_char, len=*), parameter :: asc_B3g_lbl = 'OITASC' character(kind=c_char, len=*), parameter :: asc_neq_B3g_lbl = 'ASCB3g' real(c_double), allocatable :: grid(:), mep(:), asc_Ag(:), asc_B3g(:), asc_neq_B3g(:), areas(:) integer(c_int) :: grid_size, irr_grid_size real(c_double) :: energy ! Reference values for scalar quantities integer(c_int), parameter :: ref_size = 576, ref_irr_size = 72 real(c_double), parameter :: ref_energy = -0.437960027982 if (.not. pcmsolver_is_compatible_library()) then write (error_unit, *) 'PCMSolver library not compatible!' stop end if ! Open a file for the output... inquire (file='Fortran_host.out', opened=log_open, & exist=log_exist) if (log_exist) then open (unit=output_unit, & file='Fortran_host.out', & status='unknown', & form='formatted', & access='sequential') close (unit=output_unit, status='delete') end if open (unit=output_unit, & file='Fortran_host.out', & status='new', & form='formatted', & access='sequential') rewind (output_unit) write (output_unit, *) 'Starting a PCMSolver calculation' nr_nuclei = 6_c_int allocate (charges(nr_nuclei)) allocate (coordinates(3*nr_nuclei)) ! Use C2H4 in D2h symmetry charges = (/6.0_c_double, 1.0_c_double, 1.0_c_double, & 6.0_c_double, 1.0_c_double, 1.0_c_double/) coordinates = (/0.0_c_double, 0.0_c_double, 1.257892_c_double, & 0.0_c_double, 1.745462_c_double, 2.342716_c_double, & 0.0_c_double, -1.745462_c_double, 2.342716_c_double, & 0.0_c_double, 0.0_c_double, -1.257892_c_double, & 0.0_c_double, 1.745462_c_double, -2.342716_c_double, & 0.0_c_double, -1.745462_c_double, -2.342716_c_double/) ! This means the molecular point group has three generators: ! the Oxy, Oxz and Oyz planes symmetry_info = (/3, 4, 2, 1/) host_input = pcmsolver_input() pcm_context = pcmsolver_new(PCMSOLVER_READER_HOST, & nr_nuclei, charges, coordinates, & symmetry_info, host_input, & c_funloc(host_writer)) call pcmsolver_print(pcm_context) grid_size = pcmsolver_get_cavity_size(pcm_context) irr_grid_size = pcmsolver_get_irreducible_cavity_size(pcm_context) allocate (grid(3*grid_size)) grid = 0.0_c_double call pcmsolver_get_centers(pcm_context, grid) allocate (areas(grid_size)) call pcmsolver_get_areas(pcm_context, areas) allocate (mep(grid_size)) mep = 0.0_c_double mep = nuclear_mep(nr_nuclei, charges, reshape(coordinates, (/3_c_int, nr_nuclei/)), & grid_size, reshape(grid, (/3_c_int, grid_size/))) call pcmsolver_set_surface_function(pcm_context, grid_size, mep, pcmsolver_fstring_to_carray(mep_lbl)) ! This is the Ag irreducible representation (totally symmetric) call pcmsolver_compute_asc(pcm_context, & pcmsolver_fstring_to_carray(mep_lbl), & pcmsolver_fstring_to_carray(asc_lbl), & irrep=0_c_int) allocate (asc_Ag(grid_size)) asc_Ag = 0.0_c_double call pcmsolver_get_surface_function(pcm_context, grid_size, asc_Ag, pcmsolver_fstring_to_carray(asc_lbl)) energy = pcmsolver_compute_polarization_energy(pcm_context, & pcmsolver_fstring_to_carray(mep_lbl), & pcmsolver_fstring_to_carray(asc_lbl)) write (output_unit, '(A, F20.12)') 'Polarization energy = ', energy allocate (asc_neq_B3g(grid_size)) asc_neq_B3g = 0.0_c_double ! This is the B3g irreducible representation call pcmsolver_compute_response_asc(pcm_context, & pcmsolver_fstring_to_carray(mep_lbl), & pcmsolver_fstring_to_carray(asc_neq_B3g_lbl), & irrep=3_c_int) call pcmsolver_get_surface_function(pcm_context, grid_size, asc_neq_B3g, pcmsolver_fstring_to_carray(asc_neq_B3g_lbl)) ! Equilibrium ASC in B3g symmetry. ! This is an internal check: the relevant segment of the vector ! should be the same as the one calculated using pcmsolver_compute_response_asc allocate (asc_B3g(grid_size)) asc_B3g = 0.0_c_double ! This is the B3g irreducible representation call pcmsolver_compute_asc(pcm_context, & pcmsolver_fstring_to_carray(mep_lbl), & pcmsolver_fstring_to_carray(asc_B3g_lbl), & irrep=3_c_int) call pcmsolver_get_surface_function(pcm_context, grid_size, asc_B3g, pcmsolver_fstring_to_carray(asc_B3g_lbl)) ! Check that everything calculated is OK ! Cavity size if (grid_size .ne. ref_size) then write (error_unit, *) 'Error in the cavity size, please file an issue on: https://github.com/PCMSolver/pcmsolver' stop else write (output_unit, *) 'Test on cavity size: PASSED' end if ! Irreducible cavity size if (irr_grid_size .ne. ref_irr_size) then write (error_unit, *) 'Error in the irreducible cavity size, please file an issue on: https://github.com/PCMSolver/pcmsolver' stop else write (output_unit, *) 'Test on irreducible cavity size: PASSED' end if ! Polarization energy if (.not. check_unsigned_error(energy, ref_energy, 1.0e-7_c_double)) then write (error_unit, *) 'Error in the polarization energy, please file an issue on: https://github.com/PCMSolver/pcmsolver' stop else write (output_unit, *) 'Test on polarization energy: PASSED' end if ! Surface functions call test_surface_functions(grid_size, mep, asc_Ag, asc_B3g, asc_neq_B3g, areas) call pcmsolver_save_surface_function(pcm_context, pcmsolver_fstring_to_carray(mep_lbl)) call pcmsolver_load_surface_function(pcm_context, pcmsolver_fstring_to_carray(mep_lbl)) call pcmsolver_write_timings(pcm_context) call pcmsolver_delete(pcm_context) deallocate (charges) deallocate (coordinates) deallocate (grid) deallocate (mep) deallocate (asc_Ag) deallocate (asc_B3g) deallocate (asc_neq_B3g) deallocate (areas) close (output_unit) end program pcm_fortran_host