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 <http://www.gnu.org/licenses/>.
!
! For information on the complete list of contributors to the
! PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
!

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