1! 2! PCMSolver, an API for the Polarizable Continuum Model 3! Copyright (C) 2020 Roberto Di Remigio, Luca Frediani and contributors. 4! 5! This file is part of PCMSolver. 6! 7! PCMSolver is free software: you can redistribute it and/or modify 8! it under the terms of the GNU Lesser General Public License as published by 9! the Free Software Foundation, either version 3 of the License, or 10! (at your option) any later version. 11! 12! PCMSolver is distributed in the hope that it will be useful, 13! but WITHOUT ANY WARRANTY; without even the implied warranty of 14! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! GNU Lesser General Public License for more details. 16! 17! You should have received a copy of the GNU Lesser General Public License 18! along with PCMSolver. If not, see <http://www.gnu.org/licenses/>. 19! 20! For information on the complete list of contributors to the 21! PCMSolver API, see: <http://pcmsolver.readthedocs.io/> 22! 23 24! 25! simple input reader for cavity generator 26! written by Krzysztof Mozgawa, 2010 27! 28! RDR, 280114. Put things in makecav.F inside here directly. 29! 30subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & 31 xtscor_, ytscor_, ztscor_, ar_, & 32 xsphcor_, ysphcor_, zsphcor_, rsph_, & 33 nts_, ntsirr_, nesfp_, addsph_, & 34 xe_, ye_, ze_, rin_, masses_, avgArea_, rsolv_, ret_, & 35 nr_gen_, gen1_, gen2_, gen3_, & 36 nvert_, vert_, centr_, isphe_, pedra_, len_pedra_) & 37 bind(C, name='generatecavity_cpp') 38 39use, intrinsic :: iso_c_binding 40use pedra_precision 41use pedra_symmetry, only: get_point_group, point_group 42use pedra_cavity, only: polyhedra_driver 43use strings, only: carray_to_fstring 44 45implicit none 46 47#include "pcm_pcmdef.inc" 48#include "pcm_mxcent.inc" 49#include "pcm_pcm.inc" 50 51integer(c_int) :: maxts_, maxsph_, maxvert_ 52real(c_double) :: xtscor_(maxts_), ytscor_(maxts_), ztscor_(maxts_) 53real(c_double) :: xsphcor_(maxts_), ysphcor_(maxts_), zsphcor_(maxts_), rsph_(maxts_) 54real(c_double) :: ar_(maxts_), xe_(maxts_), ye_(maxts_), ze_(maxts_), rin_(maxts_) 55real(c_double) :: masses_(maxts_) 56real(c_double) :: avgArea_, rsolv_, ret_ 57integer(c_int) :: nts_, ntsirr_, nesfp_, addsph_ 58integer(c_int) :: nr_gen_, gen1_, gen2_, gen3_ 59integer(c_int) :: nvert_(maxts_) 60real(c_double) :: vert_(maxts_ * 30), centr_(maxts_ * 30) 61integer(c_int) :: isphe_(maxts_) 62integer(c_int) :: len_pedra_ 63character(kind=c_char, len=1), intent(in) :: pedra_(len_pedra_+1) 64 65integer(c_int) :: i, j, k, offset 66integer(c_int) :: error_code 67integer(kind=regint_k) :: pedra_unit 68logical :: pedra_open, pedra_exist 69real(c_double), allocatable :: vert(:, :, :), centr(:, :, :) 70character(len=len_pedra_) :: pedra 71type(point_group) :: pgroup 72 73 74pedra = carray_to_fstring(pedra_) 75!lvpri = 121201_regint_k 76! The following INQUIRE statement returns whether the file named cavity.off is 77! connected in logical variable off_open, whether the file exists in logical 78! variable off_exist, and the unit number in integer(kind=regint_k) variable off_unit 79pedra_unit = 121201_regint_k 80inquire(file = pedra, opened = pedra_open, & 81 exist = pedra_exist) 82if (pedra_exist) then 83 open(unit = pedra_unit, & 84 file = pedra, & 85 status = 'unknown', & 86 form = 'formatted', & 87 access = 'sequential') 88 close(unit = pedra_unit, status = 'delete') 89end if 90open(unit = pedra_unit, & 91 file = pedra, & 92 status = 'new', & 93 form = 'formatted', & 94 access = 'sequential') 95rewind(pedra_unit) 96areats = avgArea_ 97icesph = 1 98iprpcm = 3 99call get_point_group(pedra_unit, pgroup, nr_gen_, gen1_, gen2_, gen3_) 100rsolv = rsolv_ 101! These parameters are fixed see one of the original GePol papers 102omega = 40.0d+00 103fro = 0.7d+00 104! ret is the minimum radius of added spheres 105ret = ret_ 106! nesfp is the number of initial spheres. 107! nesf is the total number of spheres: initial + added 108nesfp = nesfp_ 109do i = 1, nesfp 110 xe(i) = xe_(i) 111 ye(i) = ye_(i) 112 ze(i) = ze_(i) 113 rin(i) = rin_(i) 114 alpha(i) = 1.0d0 115end do 116! Workaround to silence compilation warning 117maxsph_ = mxsp 118maxvert_= mxver 119 120! Allocate space for the arrays containing the vertices and the centers 121! of the tesserae arcs 122allocate(vert(mxts, 10, 3)) 123vert = 0.0d0 124allocate(centr(mxts, 10, 3)) 125centr = 0.0d0 126 127nesf = nesfp 128 129call polyhedra_driver(pgroup, vert, centr, masses_, pedra_unit, error_code) 130 131! Common block dark magic, it will disappear one day... 132nts_ = nts 133ntsirr_ = ntsirr 134! Pass the number of added spheres back, to update the GePolCavity 135! object in the right way. 136! nesf: total number of spheres; nesfp: number of original spheres 137addsph_ = nesf - nesfp 138do i = 1, nts 139 xtscor_(i) = xtscor(i) 140 ytscor_(i) = ytscor(i) 141 ztscor_(i) = ztscor(i) 142 ar_(i) = as(i) 143 xsphcor_(i) = xe(isphe(i)) 144 ysphcor_(i) = ye(isphe(i)) 145 zsphcor_(i) = ze(isphe(i)) 146 rsph_(i) = re(isphe(i)) 147 nvert_(i) = nvert(i) 148! Fill the vert_ and centr_ arrays 149 do j = 1, nvert(i) 150 do k = 1, 3 151 offset = i + j * nts + k * nts * nvert(i) 152 vert_(offset) = vert(i, j, k) 153 centr_(offset) = centr(i, j, k) 154 end do 155 end do 156end do 157 158do i = 1, nesf 159 isphe_(i) = isphe(i) 160 xe_(i) = xe(i) 161 ye_(i) = ye(i) 162 ze_(i) = ze(i) 163 rin_(i) = re(i) 164end do 165 166! Clean-up 167deallocate(vert) 168deallocate(centr) 169 170write(pedra_unit, *) "Error code is ", error_code 171write(pedra_unit, *) '<<< Done with PEDRA Fortran code >>>' 172 173close(pedra_unit) 174 175end subroutine generatecavity_cpp 176