From 9c1feba1c149e30035821e77062bcdc28566d8c6 Mon Sep 17 00:00:00 2001 From: Xiaojie Wu Date: Mon, 13 May 2024 10:27:29 -0700 Subject: [PATCH] refactor solvent models (#139) * refactor solvent models * flake * resolve conflict * fixed bugs due to refactor * fixed typos * fixed one more type * add more unit test for SMD model * more unit test for SMD gradient and Hessian * added to_cpu * to_cpu && to_gpu for solvent * switch to legacy SMD model * pricise unit test for solvent models * remove .mod * address comments in PySCF * bug fix * move python SMD code to experiments * move SMD functionality to SMD module from pcm module --- .flake8 | 2 +- .gitignore | 1 + CONTRIBUTING.md | 4 +- benchmarks/cupy_helper/benchmark_cart2sph.py | 24 +- gpu4pyscf/lib/CMakeLists.txt | 3 +- gpu4pyscf/lib/cupy_helper.py | 2 +- gpu4pyscf/lib/cupy_helper/cart2sph.cu | 165 +- .../lib/gvhf/nr_jk_driver_int3c2e_pass1.cu | 17 +- .../lib/gvhf/nr_jk_driver_int3c2e_pass2.cu | 15 +- gpu4pyscf/lib/solvent/CMakeLists.txt | 28 +- gpu4pyscf/lib/solvent/mnsol.F | 2376 +++++++++++++++++ gpu4pyscf/lib/solvent/mnsol_interface.f90 | 18 + gpu4pyscf/lib/solvent/mnsol_mem.F | 98 + gpu4pyscf/lib/tests/test_cupy_helper.py | 16 + gpu4pyscf/lib/tests/test_to_gpu.py | 3 - gpu4pyscf/solvent/_attach_solvent.py | 240 +- gpu4pyscf/solvent/grad/pcm.py | 39 +- gpu4pyscf/solvent/grad/smd.py | 295 +- gpu4pyscf/solvent/grad/smd_experiment.py | 218 ++ gpu4pyscf/solvent/hessian/pcm.py | 27 +- gpu4pyscf/solvent/hessian/smd.py | 254 +- gpu4pyscf/solvent/hessian/smd_experiment.py | 193 ++ gpu4pyscf/solvent/pcm.py | 39 +- gpu4pyscf/solvent/smd.py | 343 +-- gpu4pyscf/solvent/smd_experiment.py | 272 ++ gpu4pyscf/solvent/tests/test_pcm.py | 85 +- gpu4pyscf/solvent/tests/test_pcm_grad.py | 109 +- gpu4pyscf/solvent/tests/test_pcm_hessian.py | 49 +- gpu4pyscf/solvent/tests/test_smd.py | 225 +- gpu4pyscf/solvent/tests/test_smd_grad.py | 230 +- gpu4pyscf/solvent/tests/test_smd_hessian.py | 236 +- 31 files changed, 4730 insertions(+), 896 deletions(-) create mode 100644 gpu4pyscf/lib/solvent/mnsol.F create mode 100644 gpu4pyscf/lib/solvent/mnsol_interface.f90 create mode 100644 gpu4pyscf/lib/solvent/mnsol_mem.F create mode 100644 gpu4pyscf/solvent/grad/smd_experiment.py create mode 100644 gpu4pyscf/solvent/hessian/smd_experiment.py create mode 100644 gpu4pyscf/solvent/smd_experiment.py diff --git a/.flake8 b/.flake8 index 572999e5..6b35cd9c 100644 --- a/.flake8 +++ b/.flake8 @@ -15,7 +15,7 @@ ignore = E701, E731, E741, E275, F401, C901, W391, W503, W504, W291, W292, W293 -exclude = test, tests, .git, __pycache__, build, dist, __init__.py .eggs, *.egg, .coveragerc +exclude = test, tests, .git, __pycache__, build, dist, __init__.py .eggs, *.egg, .coveragerc, mnsol.F max-line-length = 160 per-file-ignores = diff --git a/.gitignore b/.gitignore index 819b4e93..663ac15c 100644 --- a/.gitignore +++ b/.gitignore @@ -16,6 +16,7 @@ output/ *.swp *~ *.log +*.mod *.h5 **/qchem_input.in .coverage diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 87250ccf..fa6c4480 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -45,8 +45,8 @@ adding features to the library for PySCF functions running on GPU devices. * While examples or documentation are not mandatory, it is highly recommended to include examples of how to invoke the new module. -* CUDA compute capability 70 (sm_70) is required. Please avoid using features - that are only available on CUDA compute capability 80 or newer. The CUDA code +* CUDA compute capability 60 (sm_60) is required. Please avoid using features + that are only available on CUDA compute capability 70 or newer. The CUDA code should be compiled and run using CUDA 11 and CUDA 12 toolkits. Thank you for your contributions! diff --git a/benchmarks/cupy_helper/benchmark_cart2sph.py b/benchmarks/cupy_helper/benchmark_cart2sph.py index 620d72c3..b3d9fa2b 100644 --- a/benchmarks/cupy_helper/benchmark_cart2sph.py +++ b/benchmarks/cupy_helper/benchmark_cart2sph.py @@ -49,4 +49,26 @@ t_cutensor = perf_cutensor.gpu_times.mean() print('kernel:', t_kernel) print('cutensor:', t_cutensor) -print('memory bandwidth:',(a.nbytes+b.nbytes)/t_kernel/1024**3, 'GB/s') \ No newline at end of file +print('memory bandwidth:',(a.nbytes+b.nbytes)/t_kernel/1024**3, 'GB/s') + +print('benchmarking cart2sph when ang=5') +a = cupy.random.random([512,21*128,512]) +b = cupy.random.random([512,11*128,512]) +perf_kernel = profiler.benchmark(cart2sph, (a,1,5,b), n_repeat=20, n_warmup=3) +perf_cutensor = profiler.benchmark(cart2sph_cutensor, (a,1,5), n_repeat=20, n_warmup=3) +t_kernel = perf_kernel.gpu_times.mean() +t_cutensor = perf_cutensor.gpu_times.mean() +print('kernel:', t_kernel) +print('cutensor:', t_cutensor) +print('memory bandwidth:',(a.nbytes+b.nbytes)/t_kernel/1024**3, 'GB/s') + +print('benchmarking cart2sph when ang=6') +a = cupy.random.random([512,28*128,512]) +b = cupy.random.random([512,13*128,512]) +perf_kernel = profiler.benchmark(cart2sph, (a,1,6,b), n_repeat=20, n_warmup=3) +perf_cutensor = profiler.benchmark(cart2sph_cutensor, (a,1,6), n_repeat=20, n_warmup=3) +t_kernel = perf_kernel.gpu_times.mean() +t_cutensor = perf_cutensor.gpu_times.mean() +print('kernel:', t_kernel) +print('cutensor:', t_cutensor) +print('memory bandwidth:',(a.nbytes+b.nbytes)/t_kernel/1024**3, 'GB/s') diff --git a/gpu4pyscf/lib/CMakeLists.txt b/gpu4pyscf/lib/CMakeLists.txt index a3b482d6..12f23fbf 100644 --- a/gpu4pyscf/lib/CMakeLists.txt +++ b/gpu4pyscf/lib/CMakeLists.txt @@ -16,10 +16,9 @@ # along with this program. If not, see . cmake_minimum_required (VERSION 3.19 FATAL_ERROR) # 3.19 is required by cutlass -project (gpu4pyscf C CXX CUDA) +project (gpu4pyscf C CXX CUDA Fortran) set(CMAKE_C_STANDARD "99") - set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) diff --git a/gpu4pyscf/lib/cupy_helper.py b/gpu4pyscf/lib/cupy_helper.py index fe03587e..302b7c56 100644 --- a/gpu4pyscf/lib/cupy_helper.py +++ b/gpu4pyscf/lib/cupy_helper.py @@ -395,7 +395,7 @@ def hermi_triu(mat, hermi=1, inplace=True): def cart2sph_cutensor(t, axis=0, ang=1, out=None): ''' - transform 'axis' of a tensor from cartesian basis into spherical basis + transform 'axis' of a tensor from cartesian basis into spherical basis with cutensor ''' if(ang <= 1): if(out is not None): out[:] = t diff --git a/gpu4pyscf/lib/cupy_helper/cart2sph.cu b/gpu4pyscf/lib/cupy_helper/cart2sph.cu index 2da0c139..e9dae317 100644 --- a/gpu4pyscf/lib/cupy_helper/cart2sph.cu +++ b/gpu4pyscf/lib/cupy_helper/cart2sph.cu @@ -111,6 +111,166 @@ static void _cart2sph_ang4(double *cart, double *sph, int stride, int count){ sph[sph_offset+8*stride] = 0.625835735449176134 * (g0 + g10) - 3.755014412695056800 * g3; } +__global__ +static void _cart2sph_ang5(double *cart, double *sph, int stride, int count){ + int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= count){ + return; + } + int i = idx / stride; + int j = idx % stride; + int sph_offset = 11 * stride * i + j; + int cart_offset = 21 * stride * i + j; + double g0 = cart[cart_offset+0*stride]; + double g1 = cart[cart_offset+1*stride]; + double g2 = cart[cart_offset+2*stride]; + double g3 = cart[cart_offset+3*stride]; + double g4 = cart[cart_offset+4*stride]; + double g5 = cart[cart_offset+5*stride]; + double g6 = cart[cart_offset+6*stride]; + double g7 = cart[cart_offset+7*stride]; + double g8 = cart[cart_offset+8*stride]; + double g9 = cart[cart_offset+9*stride]; + double g10 = cart[cart_offset+10*stride]; + double g11 = cart[cart_offset+11*stride]; + double g12 = cart[cart_offset+12*stride]; + double g13 = cart[cart_offset+13*stride]; + double g14 = cart[cart_offset+14*stride]; + double g15 = cart[cart_offset+15*stride]; + double g16 = cart[cart_offset+16*stride]; + double g17 = cart[cart_offset+17*stride]; + double g18 = cart[cart_offset+18*stride]; + double g19 = cart[cart_offset+19*stride]; + double g20 = cart[cart_offset+20*stride]; + sph[sph_offset+0*stride] = 3.2819102842008507 * (g1 - 2.0 * g6) + 0.6563820568401701 * g15; + sph[sph_offset+1*stride] = 8.3026492595241645 * (g4 - g11); + sph[sph_offset+2*stride] = -1.4677148983057511 * g1 + 11.7417191864460086 * g8 + 0.4892382994352504 * (g15 - 2.0 * g6) + -3.9139063954820030 * g17; + sph[sph_offset+3*stride] = -4.7935367849733241 * (g4 + g11) + 9.5870735699466483 * g13; + sph[sph_offset+4*stride] = 0.4529466511956969 * (g1 + g15 + 2.0 * g6) + -5.4353598143483630 * (g8 + g17) + 3.6235732095655755 * g19; + sph[sph_offset+5*stride] = 1.7542548368013540 * (g2 + g16) + 3.5085096736027079 * g7 + -4.6780128981369442 * (g9 + g18) + 0.9356025796273888 * g20; + sph[sph_offset+6*stride] = 0.4529466511956969 * (g0 + g10 + 2.0 * g3) + -5.4353598143483630 * (g5 + g12) + 3.6235732095655755 * g14; + sph[sph_offset+7*stride] = -2.3967683924866621 * (g2 - g16) + 4.7935367849733241 * (g9 - g18); + sph[sph_offset+8*stride] = -0.4892382994352504 * (g0 - 2.0 * g3) + 3.9139063954820030 * g5 + 1.4677148983057511 * g10 + -11.7417191864460086 * g12; + sph[sph_offset+9*stride] = 2.0756623148810411 * (g2 + g16) + -12.4539738892862477 * g7; + sph[sph_offset+10*stride] = 0.6563820568401701 * g0 + 3.2819102842008507 * (g10 - 2.0*g3); +} + +__global__ +static void _cart2sph_ang6(double *cart, double *sph, int stride, int count){ + int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= count){ + return; + } + int i = idx / stride; + int j = idx % stride; + int sph_offset = 13 * stride * i + j; + int cart_offset = 28 * stride * i + j; + double g0 = cart[cart_offset+0*stride]; + double g1 = cart[cart_offset+1*stride]; + double g2 = cart[cart_offset+2*stride]; + double g3 = cart[cart_offset+3*stride]; + double g4 = cart[cart_offset+4*stride]; + double g5 = cart[cart_offset+5*stride]; + double g6 = cart[cart_offset+6*stride]; + double g7 = cart[cart_offset+7*stride]; + double g8 = cart[cart_offset+8*stride]; + double g9 = cart[cart_offset+9*stride]; + double g10 = cart[cart_offset+10*stride]; + double g11 = cart[cart_offset+11*stride]; + double g12 = cart[cart_offset+12*stride]; + double g13 = cart[cart_offset+13*stride]; + double g14 = cart[cart_offset+14*stride]; + double g15 = cart[cart_offset+15*stride]; + double g16 = cart[cart_offset+16*stride]; + double g17 = cart[cart_offset+17*stride]; + double g18 = cart[cart_offset+18*stride]; + double g19 = cart[cart_offset+19*stride]; + double g20 = cart[cart_offset+20*stride]; + double g21 = cart[cart_offset+21*stride]; + double g22 = cart[cart_offset+22*stride]; + double g23 = cart[cart_offset+23*stride]; + double g24 = cart[cart_offset+24*stride]; + double g25 = cart[cart_offset+25*stride]; + double g26 = cart[cart_offset+26*stride]; + double g27 = cart[cart_offset+27*stride]; + sph[sph_offset+0*stride] = 4.0991046311514863 * (g1 + g15) + -13.6636821038382887 * g6; + sph[sph_offset+1*stride] = 11.8330958111587634 * g4 + -23.6661916223175268 * g11 + 2.3666191622317525 * g22; + sph[sph_offset+2*stride] = -2.0182596029148963 * (g1 - g15) + 20.1825960291489679 * (g8 - g17); + sph[sph_offset+3*stride] = -8.2908473356343109 * g4 + -5.5272315570895412 * g11 + 22.1089262283581647 * g13 + 2.7636157785447706 * g22 + -7.3696420761193888 * g24; + sph[sph_offset+4*stride] = 0.9212052595149236 * (g1 + g15 + 2.0 * g6) + -14.7392841522387776 * (g8 + g17 - g19); + sph[sph_offset+5*stride] = 2.9131068125936568 * (g4 + g22) + 5.8262136251873136 * g11 + -11.6524272503746271 * (g13 + g24) + 4.6609709001498505 * g26; + sph[sph_offset+6*stride] = -0.3178460113381421 * (g0 + g21 + 3.0*g3 + 3.0*g10) + 5.7212282040865583 * (g5 + g23) + 11.4424564081731166 * g12 + -7.6283042721154111 * (g14 + g25) + 1.0171072362820548 * g27; + sph[sph_offset+7*stride] = 2.9131068125936568 * (g2 + g16) + 5.8262136251873136 * g7 + -11.6524272503746271 * (g9 + g18) + 4.6609709001498505 * g20; + sph[sph_offset+8*stride] = 0.4606026297574618 * (g0 - g10) + 0.4606026297574618 * (g3 - g21) + -7.3696420761193888 * (g5 - g14 - g23 + g25); + sph[sph_offset+9*stride] = -2.7636157785447706 * (g2 - 2.0 * g7) + 7.3696420761193888 * g9 + 8.2908473356343109 * g16 + -22.1089262283581647 * g18; + sph[sph_offset+10*stride] = -0.5045649007287241 * (g0 + g21) + 2.5228245036436201 * (g3 + g10) + 5.0456490072872420 * (g5 + g23) + -30.2738940437234518 * g12; + sph[sph_offset+11*stride] = 2.3666191622317525 * g2 + 11.8330958111587634 * (g16 - 2.0 * g7); + sph[sph_offset+12*stride] = 0.6831841051919144 * (g0 - g21) + -10.2477615778787161 * (g3 - g10); +} + +__global__ +static void _cart2sph_ang7(double *cart, double *sph, int stride, int count){ + int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= count){ + return; + } + int i = idx / stride; + int j = idx % stride; + int sph_offset = 15 * stride * i + j; + int cart_offset =36 * stride * i + j; + double g0 = cart[cart_offset+0*stride]; + double g1 = cart[cart_offset+1*stride]; + double g2 = cart[cart_offset+2*stride]; + double g3 = cart[cart_offset+3*stride]; + double g4 = cart[cart_offset+4*stride]; + double g5 = cart[cart_offset+5*stride]; + double g6 = cart[cart_offset+6*stride]; + double g7 = cart[cart_offset+7*stride]; + double g8 = cart[cart_offset+8*stride]; + double g9 = cart[cart_offset+9*stride]; + double g10 = cart[cart_offset+10*stride]; + double g11 = cart[cart_offset+11*stride]; + double g12 = cart[cart_offset+12*stride]; + double g13 = cart[cart_offset+13*stride]; + double g14 = cart[cart_offset+14*stride]; + double g15 = cart[cart_offset+15*stride]; + double g16 = cart[cart_offset+16*stride]; + double g17 = cart[cart_offset+17*stride]; + double g18 = cart[cart_offset+18*stride]; + double g19 = cart[cart_offset+19*stride]; + double g20 = cart[cart_offset+20*stride]; + double g21 = cart[cart_offset+21*stride]; + double g22 = cart[cart_offset+22*stride]; + double g23 = cart[cart_offset+23*stride]; + double g24 = cart[cart_offset+24*stride]; + double g25 = cart[cart_offset+25*stride]; + double g26 = cart[cart_offset+26*stride]; + double g27 = cart[cart_offset+27*stride]; + double g28 = cart[cart_offset+28*stride]; + double g29 = cart[cart_offset+29*stride]; + double g30 = cart[cart_offset+30*stride]; + double g31 = cart[cart_offset+31*stride]; + double g32 = cart[cart_offset+32*stride]; + double g33 = cart[cart_offset+33*stride]; + double g34 = cart[cart_offset+34*stride]; + double g35 = cart[cart_offset+35*stride]; + sph[sph_offset+0*stride] = 4.9501391276721742 * g1 + -24.7506956383608703 * g6 + 14.8504173830165218 * g15 + -0.7071627325245963 * g28; + sph[sph_offset+1*stride] = 15.8757639708114002 * (g4 + g22) + -52.9192132360380043 * g11; + sph[sph_offset+2*stride] = -2.5945778936013020 * (g1 - g6) + 31.1349347232156219 * g8 + 4.6702402084823440 * g15 + -62.2698694464312439 * g17 + -0.5189155787202604 * g28 + 6.2269869446431247 * g30; + sph[sph_offset+3*stride] = -12.4539738892862495 * (g4 - g22) + 41.5132462976208316 * (g13 - g24); + sph[sph_offset+4*stride] = 1.4081304047606462 * g1 + 2.3468840079344107 * g6 + -28.1626080952129243 * g8 + 0.4693768015868821 * (g15 - g28) + -18.7750720634752817 * g17 + 37.5501441269505705 * g19 + 9.3875360317376408 * g30 + -12.5167147089835229 * g32; + sph[sph_offset+5*stride] = 6.6379903866747414 * (g4 + g22) + 13.2759807733494828 * g11 + -35.4026153955986160 * (g13 + g24) + 21.2415692373591725 * g26; + sph[sph_offset+6*stride] = -0.4516580379125866 * (g1 + g28) + -1.3549741137377600 * (g6 + g15) + 10.8397929099020782 * (g8 + g30) + 21.6795858198041564 * (g17 - g19 - g32) + 5.7812228852811094 * g34; + sph[sph_offset+7*stride] = -2.3899496919201728 * (g2 + g29) + -7.1698490757605189 * (g7 + g16) + 14.3396981515210360 * (g9 + g31) + 28.6793963030420720 * g18 + -11.4717585212168292 * (g20 + g33) + 1.0925484305920790 * g35; + sph[sph_offset+8*stride] = -0.4516580379125866 * (g0 + g21) + -1.3549741137377600 * (g3 + g10) + 10.8397929099020782 * (g5 + g23) + 21.6795858198041564 * g12 + -21.6795858198041564 * (g14 + g25) + 5.7812228852811094 * g27; + sph[sph_offset+9*stride] = 3.3189951933373707 * (g2 + g7 - g16 - g29) + -17.7013076977993080 * (g9 - g31) + 10.6207846186795862 * (g20 - g33); + sph[sph_offset+10*stride] = 0.4693768015868821 * (g0 - g3) + -9.3875360317376408 * g5 + -2.3468840079344107 * g10 + 18.7750720634752817 * g12 + 12.5167147089835229 * g14 + -1.4081304047606462 * g21 + 28.1626080952129243 * g23 + -37.5501441269505705 * g25; + sph[sph_offset+11*stride] = -3.1134934723215624 * (g2 + g29) + 15.5674673616078110 * (g7 + g16) + 10.3783115744052079 * (g9 + g31) + -62.2698694464312439 * g18; + sph[sph_offset+12*stride] = -0.5189155787202604 * g0 + 4.6702402084823440 * g3 + 6.2269869446431247 * g5 + 2.5945778936013020 * (g10 - g21) + -62.2698694464312439 * g12 + 31.1349347232156219 * g23; + sph[sph_offset+13*stride] = 2.6459606618019000 * (g2 - g29) + -39.6894099270284997 * (g7 - g16); + sph[sph_offset+14*stride] = 0.7071627325245963 * g0 + -14.8504173830165218 * g3 + 24.7506956383608703 * g10 + -4.9501391276721742 * g21; +} + extern "C" { __host__ int cart2sph(cudaStream_t stream, double *cart_gto, double *sph_gto, int stride, int count, int ang) @@ -123,8 +283,11 @@ int cart2sph(cudaStream_t stream, double *cart_gto, double *sph_gto, int stride, case 2: _cart2sph_ang2 <<>> (cart_gto, sph_gto, stride, count); break; case 3: _cart2sph_ang3 <<>> (cart_gto, sph_gto, stride, count); break; case 4: _cart2sph_ang4 <<>> (cart_gto, sph_gto, stride, count); break; + case 5: _cart2sph_ang5 <<>> (cart_gto, sph_gto, stride, count); break; + case 6: _cart2sph_ang6 <<>> (cart_gto, sph_gto, stride, count); break; + case 7: _cart2sph_ang7 <<>> (cart_gto, sph_gto, stride, count); break; default: - fprintf(stderr, "Ang > 4 is not supported!"); + fprintf(stderr, "Ang > 7 is not supported!\n"); return 1; } diff --git a/gpu4pyscf/lib/gvhf/nr_jk_driver_int3c2e_pass1.cu b/gpu4pyscf/lib/gvhf/nr_jk_driver_int3c2e_pass1.cu index 7d1998a2..df7486df 100644 --- a/gpu4pyscf/lib/gvhf/nr_jk_driver_int3c2e_pass1.cu +++ b/gpu4pyscf/lib/gvhf/nr_jk_driver_int3c2e_pass1.cu @@ -52,7 +52,8 @@ static int GINTrun_tasks_int3c2e_pass1_j(JKMatrix *jk, BasisProdOffsets *offsets case 0b0000: GINTint3c2e_pass1_j_kernel0000<<>>(*envs, *jk, *offsets); break; case 0b0010: GINTint3c2e_pass1_j_kernel0010<<>>(*envs, *jk, *offsets); break; case 0b1000: GINTint3c2e_pass1_j_kernel1000<<>>(*envs, *jk, *offsets); break; - default: fprintf(stderr, "roots=1 type_ijkl %d\n", type_ijkl); + default: fprintf(stderr, "rys roots 1 type_ijkl %d\n", type_ijkl); + return 1; } break; case 2: GINTint3c2e_pass1_j_kernel<2, GSIZE2_INT3C> <<>>(*envs, *jk, *offsets); break; @@ -78,7 +79,7 @@ static int GINTrun_tasks_int3c2e_pass1_j(JKMatrix *jk, BasisProdOffsets *offsets extern "C" { __host__ int GINTbuild_j_int3c2e_pass1(BasisProdCache *bpcache, - double *dm, double *rhoj, + double *dm, double *rhoj, int nao, int naux, int n_dm, int *bins_locs_ij, int *bins_locs_kl, int ncp_ij, int ncp_kl) @@ -98,22 +99,22 @@ int GINTbuild_j_int3c2e_pass1(BasisProdCache *bpcache, jk.ao_offsets_j = 0; jk.ao_offsets_k = nao + 1; jk.ao_offsets_l = nao; - + int *bas_pairs_locs = bpcache->bas_pairs_locs; int *primitive_pairs_locs = bpcache->primitive_pairs_locs; cudaStream_t streams[MAX_STREAMS]; for (n = 0; n < MAX_STREAMS; n++){ checkCudaErrors(cudaStreamCreate(&streams[n])); } - + int *idx = (int *)malloc(sizeof(int) * TOT_NF * 3); - int *l_locs = (int *)malloc(sizeof(int) * (GPU_LMAX + 2)); + int *l_locs = (int *)malloc(sizeof(int) * (GPU_LMAX + 2)); GINTinit_index1d_xyz(idx, l_locs); checkCudaErrors(cudaMemcpyToSymbol(c_idx, idx, sizeof(int) * TOT_NF*3)); checkCudaErrors(cudaMemcpyToSymbol(c_l_locs, l_locs, sizeof(int) * (GPU_LMAX + 2))); free(idx); free(l_locs); - + for (int cp_ij_id = 0; cp_ij_id < ncp_ij; cp_ij_id++){ for (int k = 0; k < ncp_kl; k++, n++){ int n_stream = n % MAX_STREAMS; @@ -132,7 +133,7 @@ int GINTbuild_j_int3c2e_pass1(BasisProdCache *bpcache, int ntasks_kl = bins_locs_kl[k+1] - bins_locs_kl[k]; if (ntasks_kl <= 0) continue; if (ntasks_ij <= 0) continue; - + BasisProdOffsets offsets; offsets.ntasks_ij = ntasks_ij; offsets.ntasks_kl = ntasks_kl; @@ -151,7 +152,7 @@ int GINTbuild_j_int3c2e_pass1(BasisProdCache *bpcache, checkCudaErrors(cudaStreamSynchronize(streams[n])); checkCudaErrors(cudaStreamDestroy(streams[n])); } - + return 0; } diff --git a/gpu4pyscf/lib/gvhf/nr_jk_driver_int3c2e_pass2.cu b/gpu4pyscf/lib/gvhf/nr_jk_driver_int3c2e_pass2.cu index 9db78fd6..e251f9c5 100644 --- a/gpu4pyscf/lib/gvhf/nr_jk_driver_int3c2e_pass2.cu +++ b/gpu4pyscf/lib/gvhf/nr_jk_driver_int3c2e_pass2.cu @@ -51,7 +51,8 @@ static int GINTrun_tasks_int3c2e_pass2_j(JKMatrix *jk, BasisProdOffsets *offsets case 0b0000: GINTint3c2e_pass2_j_kernel0000<<>>(*envs, *jk, *offsets); break; case 0b0010: GINTint3c2e_pass2_j_kernel0010<<>>(*envs, *jk, *offsets); break; case 0b1000: GINTint3c2e_pass2_j_kernel1000<<>>(*envs, *jk, *offsets); break; - default: fprintf(stderr, "roots=1 type_ijkl %d\n", type_ijkl); + default: fprintf(stderr, "rys root 1 type_ijkl %d\n", type_ijkl); + return 1; } break; case 2: GINTint3c2e_pass2_j_kernel<2, GSIZE2_INT3C> <<>>(*envs, *jk, *offsets); break; @@ -77,7 +78,7 @@ static int GINTrun_tasks_int3c2e_pass2_j(JKMatrix *jk, BasisProdOffsets *offsets extern "C" { __host__ int GINTbuild_j_int3c2e_pass2(BasisProdCache *bpcache, - double *vj, double *rhoj, + double *vj, double *rhoj, int nao, int naux, int n_dm, int *bins_locs_ij, int *bins_locs_kl, int ncp_ij, int ncp_kl) @@ -97,16 +98,16 @@ int GINTbuild_j_int3c2e_pass2(BasisProdCache *bpcache, jk.ao_offsets_j = 0; jk.ao_offsets_k = nao + 1; jk.ao_offsets_l = nao; - + int *bas_pairs_locs = bpcache->bas_pairs_locs; int *primitive_pairs_locs = bpcache->primitive_pairs_locs; cudaStream_t streams[MAX_STREAMS]; for (n = 0; n < MAX_STREAMS; n++){ checkCudaErrors(cudaStreamCreate(&streams[n])); } - + int *idx = (int *)malloc(sizeof(int) * TOT_NF * 3); - int *l_locs = (int *)malloc(sizeof(int) * (GPU_LMAX + 2)); + int *l_locs = (int *)malloc(sizeof(int) * (GPU_LMAX + 2)); GINTinit_index1d_xyz(idx, l_locs); checkCudaErrors(cudaMemcpyToSymbol(c_idx, idx, sizeof(int) * TOT_NF*3)); checkCudaErrors(cudaMemcpyToSymbol(c_l_locs, l_locs, sizeof(int) * (GPU_LMAX + 2))); @@ -126,7 +127,7 @@ int GINTbuild_j_int3c2e_pass2(BasisProdCache *bpcache, if (envs.nrys_roots > 9) { return 2; } - + int ntasks_ij = bins_locs_ij[cp_ij_id+1] - bins_locs_ij[cp_ij_id]; int ntasks_kl = bins_locs_kl[k+1] - bins_locs_kl[k]; if (ntasks_kl <= 0) continue; @@ -149,7 +150,7 @@ int GINTbuild_j_int3c2e_pass2(BasisProdCache *bpcache, checkCudaErrors(cudaStreamSynchronize(streams[n])); checkCudaErrors(cudaStreamDestroy(streams[n])); } - + return 0; } diff --git a/gpu4pyscf/lib/solvent/CMakeLists.txt b/gpu4pyscf/lib/solvent/CMakeLists.txt index 644e3ffc..c6e9e6aa 100644 --- a/gpu4pyscf/lib/solvent/CMakeLists.txt +++ b/gpu4pyscf/lib/solvent/CMakeLists.txt @@ -1,6 +1,5 @@ -# gpu4pyscf is a plugin to use Nvidia GPU in PySCF package -# -# Copyright (C) 2022 Qiming Sun + +# Copyright 2024 The GPU4PySCF Authors. All Rights Reserved. # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -15,11 +14,28 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +#set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -arch=sm_30 -rdc=true") #set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -arch=sm_80") +#set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -arch=sm_60 --ptxas-options=-v") +#set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -arch=sm_50 --ptxas-options=-v -maxrregcount=255") +set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -std=legacy -fPIC") +set(CMAKE_VERBOSE_MAKEFILE ON) add_library(solvent SHARED - pcm.cu +mnsol_interface.f90 +mnsol_mem.F +mnsol.F +pcm.cu ) -set_target_properties(solvent PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}) -set_target_properties(solvent PROPERTIES CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") +#option(BUILD_SHARED_LIBS "build shared libraries" 1) +#option(ENABLE_STATIC "Enforce static library build" 0) +#if(ENABLE_STATIC) +# set(BUILD_SHARED_LIBS 0) +#endif() + +set_target_properties(solvent PROPERTIES + LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}) +set_target_properties(solvent PROPERTIES + CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") + diff --git a/gpu4pyscf/lib/solvent/mnsol.F b/gpu4pyscf/lib/solvent/mnsol.F new file mode 100644 index 00000000..1e0981de --- /dev/null +++ b/gpu4pyscf/lib/solvent/mnsol.F @@ -0,0 +1,2376 @@ +c copied from https://github.com/nwchemgit/nwchem/blob/master/src/solvation/mnsol.F +c MN solvation models --> +C +C THIS IS A PROGRAM TO CALCULATE THE CDS TERM FOR MINNESOTA SOLVATION MODEL SMD (MARENICH, CRAMER, TRUHLAR) +C JUN-14-2010 +C +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C +C +C |__________________:__ CDSSET +C | : | +C | : |______ CDSCALC +C | : | +C | : |______ SMD_CDS_AQ +C | : | +C | : |______ SMD_CDS_NAQ +C | : | +C | : |______ VDWRAD +C | : | +C | : |______ CDS_EG +C | : |_____ SMXCDS +C | : | |_____ IJ0 +C | : | +C | : |_____ DAREAL +C | : |_____ SCOPYMN +C | : |_____ IJ0 +C | : |_____ DSCALMN +C | : |_____ DOTMN +C | : | |_____ DDOTMN +C | : | +C | : |_____ DSWAPMN +C | : +C | : +C +C +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C +C * CDSSET +C CDS TERMS FOR MINNESOTA SOLVATION MODELS + SUBROUTINE CDSSET(ICDS,GCDS,AREACDS,NAT,C,IAN,DCDS,X, + & SOLA,SOLB,SOLC,SOLG,SOLH,SOLN) +C +C USE THE FOLLOWING VALUES OF ICDS: +C ICDS=1 FOR WATER +C ICDS=2 FOR ANY NONAQUEOUS SOLVENT +C IF ICDS=2 YOU NEED TO PROVIDE THE FOLLOWING SOLVENT DESCRIPTORS (see http://comp.chem.umn.edu/solvation/mnsddb.pdf): +C SOLA (Abraham's hydrogen-bond acidity paramerer), +C SOLB (Abraham's hydrogen-bond basicity paramerer), +C SOLC (carbon aromaticity), +C SOLG (macroscopic surface tension in cal/(mol A^2), +C SOLH (electronegative halogenicity), +C SOLN (refractive index) +C + IMPLICIT DOUBLE PRECISION (A-H,O-Z) +C +C CAVITY DISPERSION SOLVENT STRUCTURE (CDS) DRIVER + CALL CDSCALC(C,DCDS,IAN,ICDS,NAT,GCDS,AREACDS,X, + & SOLA,SOLB,SOLC,SOLG,SOLH,SOLN) +C + RETURN + END +C +C * CDSCALC + SUBROUTINE CDSCALC(C,DCDS,IAN,ICDS,NAT,GCDS,AREACDS,X, + & SOLA,SOLB,SOLC,SOLG,SOLH,SOLN) +C + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + LOGICAL LGR +C + DIMENSION C(3,*),DCDS(3,*),IAN(*),X(*) + DIMENSION SIGMA(150),HSIGMA(150) +C + GCDS=0.D0 + AREACDS=0.D0 +C +C DO WE NEED CDS GRADIENTS? LGR=.FALSE. OR LGR=.TRUE. +C LGR=.FALSE. + LGR=.TRUE. +C +C FINDING THE CDS PARAMETERS +C + DO I=1,150 + SIGMA(I)=0.D0 + HSIGMA(I)=0.D0 + ENDDO + CSSIGM=0.D0 +C + IF (ICDS.EQ.1) CALL SMD_CDS_AQ(SIGMA,HSIGMA,CSSIGM) + IF (ICDS.EQ.2) CALL SMD_CDS_NAQ(SIGMA,HSIGMA,CSSIGM, + & SOLA,SOLB,SOLC,SOLG,SOLH,SOLN) + IF (ICDS.NE.1.AND.ICDS.NE.2.AND.ICDS.NE.3) RETURN +C +C MEMORY ALLOCATION + +C +C CALL SOME EXTERNAL FUNCTION HERE TO DEFINE THE BEGINNING OF X(*) AS LXMEM !!!! +C + LXMEM = 0 ! set LXMEM = 0 for now XXX FIX ME XXXX + LRAD = LXMEM + 1 + LCDSA = LRAD + NAT + LAREA = LCDSA + NAT + LDATAR = LAREA + NAT + LSTS = LDATAR + 3 * NAT * NAT + LCOT = LSTS + NAT + LDSTS = LCOT + NAT * (NAT + 1) / 2 + LDCOTDR = LDSTS + 3 * NAT * NAT + LNC = LDCOTDR + NAT * (NAT + 1) / 2 + LDAREA = LNC + NAT + 1 + LRLIO = LDAREA + 3 * (NAT + 1) + LURLIO = LRLIO + NAT * (NAT + 1) / 2 + LLAB = LURLIO + 3 * NAT * NAT + LNCNCT = LLAB + NAT + LCONECT = LNCNCT + NAT * (2 * NAT + 1) + LCTHETA = LCONECT + NAT * NAT + LSTHETA = LCTHETA + NAT * NAT + LSIT = LSTHETA + NAT + LDCSIT = LSIT + NAT + LDJCOSN = LDCSIT + 3 * NAT + LCOSN = LDJCOSN + 3 * 3 * NAT * NAT + LDSTETA = LCOSN + 3 * NAT * NAT + LDCTETA = LDSTETA + 3 * NAT + LDCOSN = LDCTETA + 3 * NAT * NAT + LWORK = LDCOSN + 3 * 3 * NAT * NAT + LDIWORK = LWORK + 2 * NAT + 1 + LDJWORK = LDIWORK + 3 * NAT + LDKWORK = LDJWORK + 3 * NAT + LD0WORK = LDKWORK + 3 * NAT + LDWORKR = LD0WORK + 3 * NAT + LDCAODD = LDWORKR + NAT + LDSITR = LDCAODD + 3 * (NAT + 1) + LDCAPLY = LDSITR + NAT + LDICOSN = LDCAPLY + 3 * (NAT + 1) + LDCASLC = LDICOSN + 3 * 3 * NAT * NAT + LDCTETR = LDCASLC + 3 * (NAT + 1) + LDSTETR = LDCTETR + NAT + LDCOSNR = LDSTETR + NAT + LEND = LDCOSNR + 3 * NAT * NAT + INEED = LEND - LRAD +C +C YOU NEED TO CHECK IF YOU HAVE INEED WORDS OF MEMORY BY CALLING SOME EXTERNAL FUNCTION +C +C +C VAN DER WAALS RADII FOR SASA: Bondi's radii + 0.4Angs +C + SOLVRD=0.4D0 + CALL VDWRAD(X(LRAD),IAN,NAT,SOLVRD) +C +C CALCULATION OF CDS CONTRIBUTIONS TO ENERGY AND GRADIENT +C + CALL CDS_EG(CSSIGM,CDST,TAREA,NAT,LGR, + $C,DCDS,IAN,SIGMA,HSIGMA,X(LRAD), + $X(LCDSA),X(LAREA),X(LDATAR),X(LSTS),X(LCOT),X(LDSTS), + $X(LDCOTDR),X(LNC),X(LDAREA),X(LRLIO),X(LURLIO), + $X(LLAB),X(LNCNCT),X(LCONECT),X(LCTHETA),X(LSTHETA),X(LSIT), + $X(LDCSIT),X(LDJCOSN),X(LCOSN),X(LDSTETA),X(LDCTETA),X(LDCOSN), + . X(LWORK),X(LDIWORK),X(LDJWORK), + $X(LDKWORK),X(LD0WORK),X(LDWORKR),X(LDCAODD),X(LDSITR), + . X(LDCAPLY),X(LDICOSN),X(LDCASLC), + $X(LDCTETR),X(LDSTETR),X(LDCOSNR)) +C +C MEMORY RELEASE +C CALL SOME EXTERNAL FUNCTION TO RELEASE INEED WORDS OF MEMORY USED FOR X(*) +C + GCDS=CDST + AREACDS=TAREA +C +C GCDS IS A NONELECTROSTATIC CONTRIBUTION IN KCAL/MOL TO THE TOTAL ENERGY IN SOLUTION, IT MUST BE ADDED TO THE TOTAL ENERGY +C AREACDS IS THE TOTAL SURFACE AREA (FOR DEBUGGING PURPOSES ONLY) +C + RETURN + END +C +C * SMD_CDS_AQ +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C SMD_CDS_AQ CALCULATES SURFACE TENSION COEFFICIENTS +C FOR AQUEOUS SOLUTION FOR SMD +C ATOMIC SURFACE TENSION COEFFICIENTS ARE OF 3 TYPES: +C 1.) Zero order surface tension coefficients, +C SIGMA(I): atomic surface tension of atom with nuclear number I<101 +C 2.) First order modification for H, +C HSIGMA(I): H bonded to a heavy atom with nuclear number I +C 3.) Other modifications, +C SIGMA(I+100): modification for atom X bonded to Y. +C I is the atomic pair index shown as following: +C 1: C-C (1) +C 2: C-C (2) +C 3: O-C +C 4: O-O +C 5: N-C (1) +C 6: O-N +C 7: S-S +C 8: ... +C 9: ... +C 10: C-N +C 11: N-C (2) +C 12: H-(NN) +C 13: H-(OH) +C 14: O-P +C 15: S-P +C 16: N-C (triple bonded to C) +C 17: O-Si +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + SUBROUTINE SMD_CDS_AQ(SIGMA,HSIGMA,CSSIGM) +C + IMPLICIT DOUBLE PRECISION (A-H,O-Z) +C + DIMENSION SIGMA(150),HSIGMA(150) + DIMENSION SIGMA_DATA(150),HSIGMA_DATA(150) +C + DATA SIGMA_DATA/ + . 48.69, 0.00, 0.00, 0.00, 0.00, + . 129.74, 0.00, 0.00, 38.18, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . -9.10, 9.82, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, -8.72, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . -72.95, 0.00, 68.69, 0.00, -48.22, + . 121.98, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 68.85, 0.00, + . 84.10, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00/ + DATA HSIGMA_DATA/ + . 5*0.00, -60.77, 0.00, 0.00, + $ 7*0.00, 0.00,134*0.00/ +C + CSSIGM=0.D0 +C + DO I=1,150 + SIGMA(I)=SIGMA_DATA(I) + HSIGMA(I)=HSIGMA_DATA(I) + ENDDO +C + RETURN + END +C +C * SMD_CDS_NAQ +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C SMD_CDS_NAQ CALCULATES SURFACE TENSION COEFFICIENTS +C FOR NONAQUEOUS SOLUTION FOR SMD +C ATOMIC SURFACE TENSION COEFFICIENTS ARE OF 3 TYPES: +C 1.) Zero order surface tension coefficients, +C SIGMA_x(I): atomic surface tension of atom with nuclear number I<101 +C where if +C x = N is multiplied by SolN +C x = A is multiplied by SolA +C x = B is multiplied by SolB +C 2.) First order modification for H, +C HSIGMA_x(I): H bonded to a heavy atom with nuclear number I +C where if +C x = N is multiplied by SolN +C x = A is multiplied by SolA +C x = B is multiplied by SolB +C 3.) Other modifications, +C SIGMA_x(I+100): modification for atom X bonded to Y. +C I is the atomic pair index shown as following: +C 1: C-C (1) +C 2: C-C (2) +C 3: O-C +C 4: O-O +C 5: N-C (1) +C 6: O-N +C 7: S-S +C 8: ... +C 9: ... +C 10: C-N +C 11: N-C (2) +C 12: H-(NN) +C 13: H-(OH) +C 14: O-P +C 15: S-P +C 16: N-C (triple bonded to C) +C 17: O-Si +C where if +C x = N is multiplied by SolN +C x = A is multiplied by SolA +C x = B is multiplied by SolB +C MOLECULAR SURFACE TENSION COEFFICIENTS ARE STORED IN THE DATA +C ARRAY SIGMA_MOL(I), where +C element 1 is multiplied by SolG +C element 2 is multiplied by (SolB)**2 +C element 3 is multiplied by (SolC)**2 +C element 4 is multiplied by (SolH)**2 +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + SUBROUTINE SMD_CDS_NAQ(SIGMA,HSIGMA,CSSIGM, + & SOLA,SOLB,SOLC,SOLG,SOLH,SOLN) +C + IMPLICIT DOUBLE PRECISION (A-H,O-Z) +C +C YOU NEED TO PROVIDE THE FOLLOWING SOLVENT DESCRIPTORS (see http://comp.chem.umn.edu/solvation/mnsddb.pdf): +C SOLA (Abraham's hydrogen-bond acidity paramerer), +C SOLB (Abraham's hydrogen-bond basicity paramerer), +C SOLC (carbon aromaticity), +C SOLG (macroscopic surface tension in cal/(mol A^2), +C SOLH (electronegative halogenicity), +C SOLN (refractive index) + DIMENSION SIGMA(150),HSIGMA(150) + DIMENSION SIGMA_N(150),SIGMA_A(150),SIGMA_B(150) + DIMENSION HSIGMA_N(150),HSIGMA_A(150),HSIGMA_B(150) + DIMENSION SIGMA_MOL(4) + DATA SIGMA_N/ + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 58.10, 32.62, -17.56, 0.00, 0.00, + . 0.00, 0.00, 0.00, -18.04, 0.00, + . -33.17, -24.31, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, -35.42, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . -62.05, 0.00, -15.70, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, -99.76, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00/ + DATA SIGMA_A/ + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 48.10, 0.00, 193.06, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 95.99, 0.00, -41.00, + . 0.00, 0.00, 0.00, 0.00, 152.20, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00/ + DATA SIGMA_B/ + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 32.87, 0.00, -43.79, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, -128.16, 0.00, + . 79.13, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00, + . 0.00, 0.00, 0.00, 0.00, 0.00/ + DATA HSIGMA_N/ + . 5*0.00,-36.37, 0.00,-19.39,7*0.00,0.00,134*0.00/ + DATA HSIGMA_A/ + . 150*0.00/ + DATA HSIGMA_B/ + . 150*0.00/ + DATA SIGMA_MOL/ + . 0.35, 0.00, -4.19, -6.68/ +C Zero order atomic surface tensions SIGMA_x(1)-SIGMA_x(100) +C First and higher order modification SIGMA_x(101)-SIGMA_x(150) +C + DO I=1,150 + SIGMA(I)=SIGMA_N(I)*SOLN+SIGMA_A(I)*SOLA+SIGMA_B(I)*SOLB + HSIGMA(I)=HSIGMA_N(I)*SOLN+HSIGMA_A(I)*SOLA+HSIGMA_B(I)*SOLB + END DO +C + CSSIGM=SIGMA_MOL(1)*SOLG+SIGMA_MOL(2)*SOLB*SOLB+ + $ SIGMA_MOL(3)*SOLC*SOLC+SIGMA_MOL(4)*SOLH*SOLH +C + RETURN + END +C +C * VDWRAD +C VAN DER WAALS RADII FOR CDS TERM + SUBROUTINE VDWRAD(RAD,IAN,NAT,SOLVRD) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION RAD(*),IAN(*) + DIMENSION BONDI(102) +C +C SOURCE (except H): +C Van der Waals radii of Mantina et al. +C (Mantina, M.; Valero, R.; Cramer, C. J.; Truhlar, D. G. +C "Atomic Radii of the Elements." In CRC Handbook of Chemistry +C and Physics, 91st Edition, 2010-2011; Haynes, W. M., Ed.; +C CRC Press: Boca Raton, FL, 2010; pp 9-49 -- 9-50) +C +C SOURCE (H): +C Bondi's radii (Bondi, A. J. Phys. Chem. 1964, 68, 441) +C + DATA BONDI/ + $ 1.20, ! H + $ 1.40, ! He + $ 1.82, ! Li + $ 1.53, ! Be + $ 1.92, ! B + $ 1.70, ! C + $ 1.55, ! N + $ 1.52, ! O + $ 1.47, ! F + $ 1.54, ! Ne + $ 2.27, ! Na + $ 1.73, ! Mg + $ 1.84, ! Al + $ 2.10, ! Si + $ 1.80, ! P + $ 1.80, ! S + $ 1.75, ! Cl + $ 1.88, ! Ar + $ 2.75, ! K + $ 2.31, ! Ca + $ 2.16, ! Sc + $ 1.87, ! Ti + $ 1.79, ! V + $ 1.89, ! Cr + $ 1.97, ! Mn + $ 1.94, ! Fe + $ 1.92, ! Co + $ 1.84, ! Ni + $ 1.86, ! Cu + $ 2.10, ! Zn + $ 1.87, ! Ga + $ 2.11, ! Ge + $ 1.85, ! As + $ 1.90, ! Se + $ 1.85, ! Br + $ 2.02, ! Kr + $ 3.03, ! Rb + $ 2.49, ! Sr + $ 2.19, ! Y + $ 1.86, ! Zr + $ 2.07, ! Nb + $ 2.09, ! Mo + $ 2.09, ! Tc + $ 2.07, ! Ru + $ 1.95, ! Rh + $ 2.02, ! Pd + $ 2.03, ! Ag + $ 2.30, ! Cd + $ 1.93, ! In + $ 2.17, ! Sn + $ 2.06, ! Sb + $ 2.06, ! Te + $ 1.98, ! I + $ 2.16, ! Xe + $ 3.43, ! Cs + $ 2.68, ! Ba + $ 2.40, ! La + $ 2.35, ! Ce + $ 2.39, ! Pr + $ 2.29, ! Nd + $ 2.36, ! Pm + $ 2.29, ! Sm + $ 2.33, ! Eu + $ 2.37, ! Gd + $ 2.21, ! Tb + $ 2.29, ! Dy + $ 2.16, ! Ho + $ 2.35, ! Er + $ 2.27, ! Tm + $ 2.42, ! Yb + $ 2.21, ! Lu + $ 2.12, ! Hf + $ 2.17, ! Ta + $ 2.10, ! W + $ 2.17, ! Re + $ 2.16, ! Os + $ 2.02, ! Ir + $ 2.09, ! Pt + $ 2.17, ! Au + $ 2.09, ! Hg + $ 1.96, ! Tl + $ 2.02, ! Pb + $ 2.07, ! Bi + $ 1.97, ! Po + $ 2.02, ! At + $ 2.20, ! Rn + $ 3.48, ! Fr + $ 2.83, ! Ra + $ 2.60, ! Ac + $ 2.37, ! Th + $ 2.43, ! Pa + $ 2.40, ! U + $ 2.21, ! Np + $ 2.43, ! Pu + $ 2.44, ! Am + $ 2.45, ! Cm + $ 2.44, ! Bk + $ 2.45, ! Cf + $ 2.45, ! Es + $ 2.45, ! Fm + $ 2.46, ! Md + $ 2.46/ ! No +C + DO I=1,NAT + RAD(I)=0.D0 + DO K=1,102 + IF (IAN(I).EQ.K) THEN + RAD(I)=BONDI(K)+SOLVRD + ENDIF + ENDDO + ENDDO + RETURN + END +C +C * CDS_EG + SUBROUTINE CDS_EG(CSSIGM,CDST,TAREA,NAT,LGR, + $C,DCDS,IAN,SIGMA,HSIGMA,RAD, + $CDSA,AREA,DATAR,STS,COT,DSTS,DCOTDR, + $NC,DAREA,RLIO,URLIO, + $LAB,NCNCT,CONECT,CTHETA,STHETA,SIT, + $DCSIT,DJCOSN,COSN,DSTETA,DCTETA,DCOSN,WORK,DIWORK,DJWORK, + $DKWORK,D0WORK,DWORKR,DCAODD,DSITR,DCAPLY,DICOSN,DCASLC, + $DCTETR,DSTETR,DCOSNR) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (TOANGS=0.52917724924D+00,TOKCAL=627.509451D+00) + LOGICAL LGR,CONECT + DIMENSION C(3,*),DCDS(3,*),IAN(*) + DIMENSION SIGMA(150),HSIGMA(150) +C + DIMENSION CDSA(*),AREA(*) + DIMENSION DATAR(3,NAT,*) +C + DIMENSION STS(*),COT(*) + DIMENSION DSTS(3,NAT,*),DCOTDR(*) +C + DIMENSION RAD(*),NC(0:*),DAREA(3,0:*) + DIMENSION RLIO(*),URLIO(3,NAT,*) + DIMENSION LAB(*),NCNCT(2*NAT+1,*) + 1 ,CONECT(NAT,*) + DIMENSION CTHETA(NAT,*),STHETA(*),SIT(*) + 1 ,DCSIT(3,*),DJCOSN(3,3,NAT,*) + 2 ,COSN(3,NAT,*),DSTETA(3,*) + 3 ,DCTETA(3,NAT,*),DCOSN(3,3,NAT,*) + 4 ,WORK(*),DIWORK(3,*),DJWORK(3,*) + 5 ,DKWORK(3,*),D0WORK(3,*),DWORKR(*) + 6 ,DCAODD(3,0:*),DSITR(*) + 6 ,DCAPLY(3,0:*) + 7 ,DICOSN(3,3,NAT,*),DCASLC(3,0:*) + 8 ,DCTETR(*),DSTETR(*) + 9 ,DCOSNR(3,NAT,*) +C + DO I=1,NAT + CDSA(I)=0.D0 + AREA(I)=0.D0 + STS(I)=0.D0 + DO J=1,NAT + DO K=1,3 + DATAR(K,I,J)=0.D0 + ENDDO + ENDDO + ENDDO +C +C MATRIX OF INTERATOMIC DISTANCES IN ANGSTROMS +C + K12=0 + DO K1=1,NAT + DO K2=1,K1 + K12=K12+1 + IF (K1.EQ.K2) THEN + RLIO(K12)=0.D0 + URLIO(1,K1,K2)=0.D0 + URLIO(2,K1,K2)=0.D0 + URLIO(3,K1,K2)=0.D0 + ELSE + X=(C(1,K1)-C(1,K2))*TOANGS + Y=(C(2,K1)-C(2,K2))*TOANGS + Z=(C(3,K1)-C(3,K2))*TOANGS + R=DSQRT(X*X+Y*Y+Z*Z) + URLIO(1,K1,K2)=-X/R + URLIO(2,K1,K2)=-Y/R + URLIO(3,K1,K2)=-Z/R + URLIO(1,K2,K1)=X/R + URLIO(2,K2,K1)=Y/R + URLIO(3,K2,K1)=Z/R + RLIO(K12)=R + ENDIF + ENDDO + ENDDO +C +C CALCULATION OF CDS CONTRIBUTIONS TO ENERGY +C + CALL SMXCDS + $(IAN,SIGMA,HSIGMA,STS,COT,DSTS,DCOTDR,RLIO,URLIO, + $NAT,LGR) +C + TAREA=0.D0 + CDST=0.D0 + DO K=1,NAT + AREA0=0.D0 + NCROSS=0 + CALL DAREAL + $(NAT,AREA0,NCROSS,K,LGR, + $RAD,NC,DAREA,RLIO,URLIO,LAB,NCNCT,CONECT,CTHETA,STHETA,SIT, + $DCSIT,DJCOSN,COSN,DSTETA,DCTETA,DCOSN,WORK,DIWORK,DJWORK, + $DKWORK,D0WORK,DWORKR,DCAODD,DSITR,DCAPLY,DICOSN,DCASLC, + $DCTETR,DSTETR,DCOSNR) + RAS2=RAD(K)*RAD(K) + AREA(K)=AREA0*RAS2 + CDST=CDST+AREA(K)*(STS(K)+CSSIGM)*0.001D0 + TAREA=TAREA+AREA(K) + CDSA(K)=AREA(K)*(STS(K)+CSSIGM)*0.001D0 + IF (LGR) THEN + DO L=0,NCROSS + DATAR(1,NC(L),K)=DAREA(1,L)*RAS2 + DATAR(2,NC(L),K)=DAREA(2,L)*RAS2 + DATAR(3,NC(L),K)=DAREA(3,L)*RAS2 + ENDDO + ENDIF + ENDDO +C ----- COMPUTATION OF CDS GRADIENTS ----- +C + IF (LGR) THEN + DO IAT=1,NAT + DCDSX=0.D0 + DCDSY=0.D0 + DCDSZ=0.D0 + DO J=1,NAT + DCDSX=DCDSX + $ +(STS(J)+CSSIGM)*DATAR(1,IAT,J)*0.001D0 + $ +DSTS(1,IAT,J)*AREA(J)*0.001D0 + DCDSY=DCDSY + $ +(STS(J)+CSSIGM)*DATAR(2,IAT,J)*0.001D0 + $ +DSTS(2,IAT,J)*AREA(J)*0.001D0 + DCDSZ=DCDSZ + $ +(STS(J)+CSSIGM)*DATAR(3,IAT,J)*0.001D0 + $ +DSTS(3,IAT,J)*AREA(J)*0.001D0 + ENDDO + DCDSX=DCDSX/TOKCAL*TOANGS + DCDSY=DCDSY/TOKCAL*TOANGS + DCDSZ=DCDSZ/TOKCAL*TOANGS + DCDS(1,IAT)=DCDSX + DCDS(2,IAT)=DCDSY + DCDS(3,IAT)=DCDSZ + ENDDO + ENDIF +C + RETURN + END +C +C * SMXCDS + SUBROUTINE SMXCDS + $(IAN,SIGMA,HSIGMA,STS,COT,DSTS,DCOTDR,RLIO,URLIO, + $NAT,LGR) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + LOGICAL LGR + DIMENSION IAN(*) + DIMENSION SIGMA(150),HSIGMA(150) + DIMENSION NATCNV(100),RKKVAL(15,15) + DIMENSION STS(*),COT(*) + DIMENSION DSTS(3,NAT,*),DCOTDR(*), + $RLIO(*),URLIO(3,NAT,*) +C H C N O F Si P S Cl Br I + DATA NATCNV/1,4*0,2,3,4,5,4*0,11,9,6,7,17*0,8,17*0,10,47*0/ +C +C --- RKKVAL EXPANDED TO 15 X 15 DATA ARRAY --- +C +C --- TYPE X REFERS TO : --- +C +C Hydrogen - type 1 +C Carbon - type 2 +C Nitrogen - type 3 +C Oxygen - type 4 +C Fluorine - type 5 +C Sulfur - type 6 +C Chlorine - type 7 +C Bromine - type 8 +C Phos. - type 9 +C Iodine - type 10 +C Silicon - type 11 +C +C + DATA ((RKKVAL(J,I), I=1,15), J=1,15) +C +C --- J => TYPES DEFINED ABOVE --- +C --- I => NATCNV ARRAY ELEMENTS --- +C +C Set for type 1 + . /0.00D0, 1.55D0, 1.55D0, 1.55D0, 0.00D0, + . 2.14D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 2 + . 1.55D0, 1.84D0, 1.84D0, 1.84D0, 1.84D0, + . 2.20D0, 2.10D0, 2.30D0, 2.20D0, 2.60D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 3 + . 1.55D0, 1.84D0, 1.85D0, 1.50D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 4 + . 1.55D0, 1.84D0, 1.50D0, 2.75D0, 0.00D0, + . 1.71D0, 0.00D0, 0.00D0, 2.10D0, 0.00D0, + . 2.10D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 5 + . 0.00D0, 1.84D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 6 + . 2.14D0, 2.20D0, 0.00D0, 1.71D0, 0.00D0, + . 2.75D0, 0.00D0, 0.00D0, 2.50D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 7 + . 0.00D0, 2.10D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 8 + . 0.00D0, 2.30D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 9 + . 0.00D0, 2.20D0, 0.00D0, 2.10D0, 0.00D0, + . 2.50D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 10 + . 0.00D0, 2.60D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 11 + . 0.00D0, 0.00D0, 0.00D0, 2.10D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 12 + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 13 + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 14 + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, +C Set for type 15 + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0, + . 0.00D0, 0.00D0, 0.00D0, 0.00D0, 0.00D0/ +C + ISTS=6 +C + DO I=1, NAT*(NAT+1)/2 + DCOTDR(I)=0.D0 + COT(I)=0.D0 + ENDDO + DO I1=1,NAT + DO I2=1,NAT + DO I3=1,3 + DSTS(I3,I1,I2)=0.D0 + ENDDO + ENDDO + ENDDO +C +C-----COMPUTE COTS AND DERIVATIVES +C + RHLD=1.4d0 + DELTAR=0.2d0 + BCCC=1.0d0 + EXPVAL=DEXP(1.00D0) +C + DO 100 I=1,NAT + II=IAN(I) + ITPC=NATCNV(II) + DO 100 J=1,NAT + IF (I.NE.J) THEN + JJ=IAN(J) + JTPC=NATCNV(JJ) + IF(ITPC.EQ.0.OR.JTPC.EQ.0) THEN + COT(IJ0(I,J))=0.0D0 + ELSE + COT(IJ0(I,J))=0.0D0 + DELTAR=0.30D0 + RHLD=RKKVAL(ITPC,JTPC) + IF((II.EQ.8.and.JJ.EQ.6).OR.(II.EQ.6.and.JJ.EQ.8)) THEN + DELTAR=0.10D0 + RHLD=1.330D0 + END IF + IF (ISTS.EQ.6) THEN + IF(II.EQ.8.and.JJ.EQ.8) THEN + DELTAR=0.30D0 + RHLD=1.80D0 + END IF + END IF + CUTOF1=RHLD+DELTAR + IF(RLIO(IJ0(I,J)).LT.CUTOF1) THEN + EXPONT=DELTAR/(RLIO(IJ0(I,J))-CUTOF1) + COT(IJ0(I,J))=DEXP(EXPONT) +C + if (LGR) + $DCOTDR(IJ0(I,J))=-COT(IJ0(I,J))*EXPONT/(RLIO(IJ0(I,J))-CUTOF1) + END IF + END IF + ENDIF + 100 CONTINUE +C +C-----CDS SIGMAS CAN BE ANY OF THESE FORMS: a, a-b, a-b(2), a-b(3) (a AND b - ARE ELEMENTS) +C-----SOME OF THE EXPRESSIONS BELOW ARE NOT USED BY THE SMD MODEL BUT KEPT HERE FOR DEBUGGING PURPOSES +C + DO 950 I=1,NAT +C +C-----ADD CDS FOR SIGMA a + STS(I)=SIGMA(IAN(I)) +C + IF (IAN(I).EQ.1) THEN +C +C-----ADD CDS FOR SIGMA H-a (USED BY SMD ONLY FOR H-C, H-O) + DO 150 J=1,NAT + NTP=IAN(J) + STS(I)=STS(I)+COT(IJ0(I,J))*HSIGMA(NTP) +C-----ADD DCDS FOR SIGMA H-a + if (LGR) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +HSIGMA(NTP)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -HSIGMA(NTP)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + end do + end if +C +C-----ADD CDS FOR SIGMA H-N(2); FOR H-N-N BOND (NOT USED BY SMD) + IF (NTP.EQ.7) THEN + DO 262 K=1,NAT + KTP=IAN(K) + IF (J.NE.K.AND.KTP.EQ.7) THEN + STS(I)=STS(I)+COT(IJ0(I,J))*COT(IJ0(J,K))*SIGMA(112) +C-----ADD DCDS FOR SIGMA H-N(2); FOR H-N-N BOND (NOT USED BY SMD) + if (LGR) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(112)*COT(ij0(J,K)) + . *DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(112)*COT(ij0(J,K)) + . *DCOTDR(ij0(I,J))*URLIO(L,J,I) + . +SIGMA(112)*COT(ij0(I,J)) + . *DCOTDR(ij0(J,K))*URLIO(L,K,J) + DSTS(L,K,I)=DSTS(L,K,I) + . -SIGMA(112)*COT(ij0(I,J)) + . *DCOTDR(ij0(J,K))*URLIO(L,K,J) + end do + end if + ENDIF + 262 CONTINUE +C +C-----ADD CDS FOR SIGMA H-O(2); FOR H-O-H BOND (NOT USED BY SMD) + ELSE IF (NTP.EQ.8) THEN + DO 263 K=1,NAT + KTP=IAN(K) + IF (I.NE.K.AND.KTP.EQ.1) THEN + STS(I)=STS(I)+COT(IJ0(I,J))*COT(IJ0(J,K))*SIGMA(113) +C-----ADD DCDS FOR SIGMA H-O(2); FOR H-O-H BOND (NOT USED BY SMD) + if (LGR) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(113)*COT(ij0(J,K)) + . *DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(113)*COT(ij0(J,K)) + . *DCOTDR(ij0(I,J))*URLIO(L,J,I) + . +SIGMA(113)*COT(ij0(I,J)) + . *DCOTDR(ij0(J,K))*URLIO(L,K,J) + DSTS(L,K,I)=DSTS(L,K,I) + . -SIGMA(113)*COT(ij0(I,J)) + . *DCOTDR(ij0(J,K))*URLIO(L,K,J) + end do + end if + ENDIF + 263 CONTINUE + ENDIF + 150 CONTINUE + ENDIF +C + IF (IAN(I).EQ.8) THEN + RTKK=0.0D0 + KOO=0 + DO 200 J=1,NAT + NTP=IAN(J) +C +C-----ADD CDS FOR SIGMA O-C + IF(NTP.EQ.6)THEN + STS(I)=STS(I)+COT(IJ0(I,J))*SIGMA(103) +C-----ADD DCDS FOR SIGMA O-C + if (LGR) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(103)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(103)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + end do + end if +C +C-----ADD CDS FOR SIGMA O-N + ELSE IF(NTP.EQ.7) THEN + STS(I)=STS(I)+COT(IJ0(I,J))*SIGMA(106) +C-----ADD DCDS FOR SIGMA O-N + if (LGR) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(106)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(106)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + end do + end if +C +C-----ADD CDS FOR SIGMA O-S (NOT USED BY SMD) + ELSE IF(NTP.EQ.16) THEN + STS(I)=STS(I)+COT(IJ0(I,J))*SIGMA(118) +C-----ADD DCDS FOR SIGMA O-S (NOT USED BY SMD) + if (LGR) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(118)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(118)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + end do + end if +C +C-----ADD CDS FOR SIGMA O-SI (NOT USED BY SMD) + Else If (NTP .EQ. 14) Then + STS(I)=STS(I)+COT(IJ0(I,J))*SIGMA(117) +C-----ADD DCDS FOR SIGMA O-SI (NOT USED BY SMD) + If (LGR) Then + Do L = 1,3 + DSTS(L,I,I)=DSTS(L,I,I) + & +SIGMA(117)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + & -SIGMA(117)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + End Do + End If +C +C-----ADD CDS FOR SIGMA O-P + ELSE IF (NTP.EQ.15) THEN + STS(I)=STS(I)+COT(IJ0(I,J))*SIGMA(114) +C-----ADD DCDS FOR SIGMA O-P + if (LGR) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(114)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(114)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + end do + end if +C +C-----ADD CDS FOR SIGMA O-O (SM6 FORM ONLY) + ELSE IF(NTP.EQ.8.AND.I.NE.J) THEN + If (ISTS.eq.6) then + STS(I)=STS(I)+COT(IJ0(I,J))*SIGMA(104) +C-----ADD DCDS FOR SIGMA O-O (SM6 FORM ONLY) + if (LGR) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(104)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(104)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + end do + end if + End if + ENDIF + 200 CONTINUE + ENDIF +C +C-----ADD CDS FOR SIGMA N-C FOR -N-C- BOND AND N-C(2) FOR N-C=O BOND; N-C(2) IS NOT USED BY SMD + IF (IAN(I).EQ.7) THEN + RTKKS=0.0D0 + RTKKS2=0.D0 + DO 250 J=1,NAT + NTP=IAN(J) + IF(NTP.EQ.6)THEN + RTKK3=0.0D0 + RTKK5=0.0D0 + DO 251 K=1,NAT + IF (K.NE.I.AND.K.NE.J) THEN + KTP=IAN(K) + NKTPC=NATCNV(IAN(K)) + NTPC =NATCNV(IAN(J)) + RHLD2=RKKVAL(NKTPC,NTPC) + CUTOF2=RHLD2+0.3D0 + IF(RLIO(IJ0(J,K)).LT.CUTOF2)THEN + RTKK2=DEXP(1.0D0-(BCCC/(1.0D0- + . ((RLIO(IJ0(J,K))-RHLD2)/0.3D0))))/(EXPVAL) + IF(KTP.EQ.8) RTKK5=RTKK5+RTKK2 + RTKK3=RTKK3+RTKK2 + ENDIF + ENDIF + 251 CONTINUE + RTKKS=RTKKS+COT(IJ0(I,J))*RTKK3**2 + RTKKS2=RTKKS2+COT(IJ0(I,J))*RTKK5 + ENDIF + 250 CONTINUE + DHOLDER=(RTKKS**1.3D0) + STS(I)=STS(I)+DHOLDER*SIGMA(105) + STS(I)=STS(I)+RTKKS2*SIGMA(111) +C-----ADD DCDS FOR SIGMA N-C FOR -N-C- BOND AND N-C(2) FOR N-C=O BOND; N-C(2) IS NOT USED BY SMD + if (LGR) then + do J=1,NAT + NTP=IAN(J) + if (NTP.EQ.6) then + SCOTC=0.0D0 + SCOTO=0.0D0 + do K=1,NAT + if ((K.ne.I).and.(K.ne.J)) then + KTP=IAN(K) + COTJK=0.0D0 + DCOTJK=0.0D0 + NKTPC=NATCNV(IAN(K)) + NTPC =NATCNV(IAN(J)) + RHLD2=RKKVAL(NKTPC,NTPC) + CUTOF2=RHLD2+0.30 + if (RLIO(ij0(J,K)).lt.CUTOF2) then + COTJK=exp(1.0D0-(BCCC/(1.0D0 + . -((RLIO(ij0(J,K))-RHLD2)/0.30))))/(EXPVAL) + DCOTJK=-COTJK*0.3/(RLIO(ij0(J,K))-CUTOF2)**2 + end if + if (KTP.eq.8) then + SCOTO=SCOTO+COTJK + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(111)*COTJK + . *DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(111)*COTJK + . *DCOTDR(ij0(I,J))*URLIO(L,J,I) + . +SIGMA(111)*COT(ij0(I,J)) + . *DCOTJK *URLIO(L,K,J) + DSTS(L,K,I)=DSTS(L,K,I) + . -SIGMA(111)*COT(ij0(I,J)) + . *DCOTJK *URLIO(L,K,J) + end do + end if + SCOTC=SCOTC+COTJK + end if + end do +C + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(105)*RTKKS**0.3D0*1.3D0 + . *SCOTC**2*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(105)*RTKKS**0.3D0*1.3D0 + . *SCOTC**2*DCOTDR(ij0(I,J))*URLIO(L,J,I) + end do +C + do K=1,NAT + if ((K.ne.I).and.(K.ne.J)) then + KTP=IAN(K) + COTJK=0.0D0 + DCOTJK=0.0D0 + NKTPC=NATCNV(IAN(K)) + NTPC =NATCNV(IAN(J)) + RHLD2=RKKVAL(NKTPC,NTPC) + CUTOF2=RHLD2+0.30 + if (RLIO(ij0(J,K)).lt.CUTOF2) then + COTJK=exp(1.0D0-(1.0D0/(1.0D0 + . -((RLIO(ij0(J,K))-RHLD2)/0.30))))/(EXPVAL) + DCOTJK=-COTJK*0.3/(RLIO(ij0(J,K))-CUTOF2)**2 + end if + do L=1,3 + DSTS(L,J,I)=DSTS(L,J,I)+SIGMA(105) + . *RTKKS**0.3D0*1.3D0 + . *2.0D0*SCOTC*COT(ij0(I,J)) + . *DCOTJK *URLIO(L,K,J) + DSTS(L,K,I)=DSTS(L,K,I)-SIGMA(105) + . *RTKKS**0.3D0*1.3D0 + . *2.0D0*SCOTC*COT(ij0(I,J)) + . *DCOTJK *URLIO(L,K,J) + end do + end if + end do + end if + end do + end if +C +C-----ADD CDS FOR SIGMA N-C(3) FOR THE N-C TRIPLE BOND + DRTKKU=0.D0 + RTKK=0.0D0 + DELTAR=0.065D0 + DO 461 J=1,NAT + NTP=IAN(J) + NTPC=NATCNV(NTP) + IF(J.EQ.I) GOTO 461 + IF(NTP.EQ.6) THEN + RHLD=1.225D0 + CUTOF1=RHLD+DELTAR + IF(RLIO(IJ0(I,J)).LT.CUTOF1) THEN + RTKK=RTKK+DEXP(1.0D0-(BCCC/(1.0D0-((RLIO(IJ0(I,J))- + . RHLD)/DELTAR))))/(EXPVAL) +C-----ADD DCDS FOR SIGMA N-C(3) FOR THE N-C TRIPLE BOND + if (LGR) then + EXPONT=DELTAR/(RLIO(ij0(I,J))-CUTOF1) + DRTKK=-DEXP(EXPONT)*EXPONT/(RLIO(ij0(I,J))-CUTOF1) + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(116)*DRTKK*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(116)*DRTKK*URLIO(L,J,I) + end do + end if + END IF + END IF + 461 CONTINUE + STS(I)=STS(I)+RTKK*SIGMA(116) + ENDIF +C +C-----ADD CDS FOR SIGMA C-C AND C-C(2); C-C(2) IS NOT USED BY SMD + IF (IAN(I).EQ.6) THEN + ITPC=NATCNV(6) + DO 500 KCCC=1,2 + RTKK=0.0D0 + DO 450 J=1,NAT + IF(J.EQ.I) GOTO 450 + NTP=IAN(J) + IF(NTP.EQ.6)THEN + NTPC=NATCNV(NTP) + IF (KCCC.EQ.1) THEN + RHLD=RKKVAL(ITPC,NTPC) + DELTAR=0.30D0 + ELSE + RHLD=1.27D0 + DELTAR=0.07D0 + ENDIF + CUTOF1=RHLD+DELTAR + IF(RLIO(IJ0(I,J)).LT.CUTOF1) THEN + RTKK=RTKK+DEXP(1.0D0-(BCCC/(1.0D0-((RLIO(IJ0(I,J))- + . RHLD)/DELTAR))))/(EXPVAL) +C-----ADD DCDS FOR SIGMA C-C AND C-C(2); C-C(2) IS NOT USED BY SMD + if (LGR) then + EXPONT=DELTAR/(RLIO(ij0(I,J))-CUTOF1) + DRTKK=-DEXP(EXPONT)*EXPONT/(RLIO(ij0(I,J))-CUTOF1) + if(KCCC.EQ.1)then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(101)*DRTKK*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(101)*DRTKK*URLIO(L,J,I) + end do + else if (KCCC.EQ.2) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(102)*DRTKK*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(102)*DRTKK*URLIO(L,J,I) + end do + end if + end if + ENDIF + ENDIF + 450 CONTINUE + IF(KCCC.EQ.1)THEN + STS(I)=STS(I)+RTKK*SIGMA(101) + ELSE IF (KCCC.EQ.2) THEN + STS(I)=STS(I)+RTKK*SIGMA(102) + ENDIF + 500 CONTINUE +C +C-----ADD CDS FOR SIGMA C-N FOR ANY C-N BOND + RTKK=0.0D0 + DO 460 J=1,NAT + NTP=IAN(J) + IF(NTP.EQ.7)THEN + RTKK=RTKK+COT(IJ0(I,J)) + ENDIF + 460 CONTINUE + HOLDRK=RTKK*RTKK + STS(I)=STS(I)+HOLDRK*SIGMA(110) +C-----ADD DCDS FOR SIGMA C-N FOR ANY C-N BOND + if (LGR) then + do J=1,NAT + NTP=IAN(J) + if (NTP.eq.7) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(110)*2.0D0 + . *RTKK*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(110)*2.0D0 + . *RTKK*DCOTDR(ij0(I,J))*URLIO(L,J,I) + end do + end if + end do + end if +C + ENDIF +C +C-----ADD CDS FOR SIGMA S-S FOR S-S BOND (NOT USED BY SMD) + IF (IAN(I).EQ.16) THEN + DO 550 J=1,NAT + NTP=IAN(J) + IF(NTP.EQ.16.AND.J.NE.I) THEN + STS(I)=STS(I)+COT(IJ0(I,J))*SIGMA(107) +C-----ADD DCDS FOR SIGMA S-S FOR S-S BOND (NOT USED BY SMD) + if (LGR) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(107)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(107)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + end do + end if +C + ENDIF + 550 CONTINUE + ENDIF +C +C-----ADD CDS FOR SIGMA S-P FOR S-P BOND (NOT USED BY SMD) + IF (IAN(I).EQ.16) THEN + DO 220 J=1,NAT + NTP=IAN(J) + IF(NTP.EQ.15)THEN + STS(I)=STS(I)+COT(IJ0(I,J))*SIGMA(115) +C-----ADD DCDS FOR SIGMA S-P FOR S-P BOND (NOT USED BY SMD) + if (LGR) then + do L=1,3 + DSTS(L,I,I)=DSTS(L,I,I) + . +SIGMA(115)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + DSTS(L,J,I)=DSTS(L,J,I) + . -SIGMA(115)*DCOTDR(ij0(I,J))*URLIO(L,J,I) + end do + end if +C + ENDIF + 220 CONTINUE + ENDIF + 950 CONTINUE + RETURN + END +C +C * IJ0 + FUNCTION IJ0(I,J) + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + IF(I.GT.J) THEN + IJ0=I*(I-1)/2+J + ELSE + IJ0=J*(J-1)/2+I + END IF + RETURN + END +C +C * DAREAL +C ******************** DAREAL SUBROUTINE ********************** +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C +C ACCESSIBLE SOLID ANGLE OF SPHERES (RADIANS). ANALYTICAL. +C ALLOW MULTIPLE CONNECTED COMPONENTS ON A SPHERE. +C REGULARIZE RARE EVENTS (POINTS SHARED BY MORE THAN 3 SPHERES) +C BY A SMALL PERTURBATION ON RADIUS. +C +C INPUT +C RAD(I) : RADII. +C NAT : NUMBER OF SPHERES. +C K : NO OF THE SPHERE STUDIED. +C LGRX : .TRUE. TO CALCULATE CARTESIAN DERIVATIVES. +C LGRR : .TRUE. TO CALCULATE DERIVATIVE VS RAD(K). +C (REF) +C RLIO(I): INTERCENTRE DISTANCES (LOW TRIANGLE), +C URLIO(3,I,J) : ASSOCIATED UNIT VECTORS: FROM i TO j. +C +C OUTPUT +C AREA : ACCESSIBLE SOLID ANGLE FOR SPHERE NO K. +C NCROSS : NUMBER OF SPHERES OVERLAPPING K. +C NC : INDICE OF SPHERES CONNECTED TO K. +C AND IF LGRX = .TRUE. : +C DAREA : DERIVATIVES OF THE ACCESSIBLE SOLID ANGLE WITH +C RESPECT TO THE CARTESIAN COORDINATES. +C THAT IS: DAREA(I,J) = dAREA(K)/dCOORD(I,NC(J)) ... J=0,NCROSS. +C AND IF LGRR = .TRUE. : +C DAREAR : DERIVATIVES OF THE ACCESSIBLE SOLID ANGLE WITH +C RESPECT TO RAD(K). +C +C By D. Liotard, December 1992. +C Derivatives by D. Rinaldi and D. Liotard +C +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + SUBROUTINE DAREAL + $(NAT,AREA,NCROSS,K,LGRX, + $RAD,NC,DAREA,RLIO,URLIO,LAB,NCNCT,CONECT,CTHETA,STHETA,SIT, + $DCSIT,DJCOSN,COSN,DSTETA,DCTETA,DCOSN,WORK,DIWORK,DJWORK, + $DKWORK,D0WORK,DWORKR,DCAODD,DSITR,DCAPLY,DICOSN,DCASLC, + $DCTETR,DSTETR,DCOSNR) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (PI=3.14159265358979D0) + LOGICAL FREEIJ,FREEJI,CONECT,LPOLY,LGRX,LGRR + DIMENSION RAD(*),NC(0:*),DAREA(3,0:*) + DIMENSION RLIO(*),URLIO(3,NAT,*) + DIMENSION LAB(*),NCNCT(2*NAT+1,*) + 1 ,CONECT(NAT,*) + DIMENSION CTHETA(NAT,*),STHETA(*),SIT(*) + 1 ,DCSIT(3,*),DJCOSN(3,3,NAT,*) + 2 ,COSN(3,NAT,*),DSTETA(3,*) + 3 ,DCTETA(3,NAT,*),DCOSN(3,3,NAT,*) + 4 ,WORK(*),DIWORK(3,*),DJWORK(3,*) + 5 ,DKWORK(3,*),D0WORK(3,*),DWORKR(*) + 6 ,DCAODD(3,0:*),DSITR(*) + 6 ,DCAPLY(3,0:*) + 7 ,DICOSN(3,3,NAT,*),DCASLC(3,0:*) + 8 ,DCTETR(*),DSTETR(*) + 9 ,DCOSNR(3,NAT,*) + DIMENSION DCNIJR(3),DVIJR(3) + A ,DCNIKR(3),DCC2I(3),DKCNIK(3,3),DAPHI(3),DBPHI(3) + B ,VIJ(3),DICOS(3,3),DJCOS(3,3),VN(3),DIVN(3,3) + C ,DJVN(3,3),DIVIJ(3,3),DJVIJ(3,3),DIWIJ(3,3) + D ,DJWIJ(3,3),DIAIJ(3),DJAIJ(3),DIBIJ(3),DJBIJ(3) + E ,DICIJ(3),DJCIJ(3),CNIJ(3),DICNIJ(3,3),DJCNIJ(3,3) + F ,CNIK(3),DICNIK(3,3) + DATA EPSI /1.D-11/ +C + MXSS=2*NAT+1 +C + PCMIN=1000 + JPCNT=0 +C + LGRR=.FALSE. +C + DASLIR=0.D0 + DCIJ=0.D0 + DSN2IJ=0.D0 + DAIJR=0.D0 + DBIJR=0.D0 + DCIJR=0.D0 + DAPOLR=0.D0 + DC2IR=0.D0 + DAODDR=0.D0 + DPHIR=0.D0 +C + DO I=1,NAT + LAB(I)=0 + STHETA(I)=0.D0 + SIT(I)=0.D0 + DWORKR(I)=0.D0 + DSITR(I)=0.D0 + DCTETR(I)=0.D0 + DSTETR(I)=0.D0 + DO J=1,3 + DCSIT(J,I)=0.D0 + DSTETA(J,I)=0.D0 + DIWORK(J,I)=0.D0 + DJWORK(J,I)=0.D0 + DKWORK(J,I)=0.D0 + D0WORK(J,I)=0.D0 + ENDDO + DO J=1,MXSS + NCNCT(J,I)=0 + ENDDO + DO J=1,NAT + CONECT(I,J)=.FALSE. + CTHETA(I,J)=0.D0 + DO KK=1,3 + COSN(KK,I,J)=0.D0 + DCTETA(KK,I,J)=0.D0 + DCOSNR(KK,I,J)=0.D0 + DO L=1,3 + DJCOSN(KK,L,I,J)=0.D0 + DCOSN(KK,L,I,J)=0.D0 + DICOSN(KK,L,I,J)=0.D0 + ENDDO + ENDDO + ENDDO + ENDDO + DO I=1,MXSS + WORK(I)=0.D0 + ENDDO + DO I=0,NAT + DO J=1,3 + DCAODD(J,I)=0.D0 + DCAPLY(J,I)=0.D0 + DCASLC(J,I)=0.D0 + ENDDO + ENDDO +C +C EPSI:THRESHOLD TO AVOID DRAMATIC CONSEQUENCES OF ROUND-OFF ERRORS, +C MUST BE GREATER THAN THE POOREST OF RELATIVE ERRORS ON FUNCTIONS: +C SQRT, ATAN. +C + TWOPI=PI+PI + FOURPI=TWOPI+TWOPI + RK=RAD(K) + NC(0) = K + NCROSS = 0 + IF (LGRX) CALL SCOPYMN(3,0D0,0,DAREA(1,0),1) +C + IF (LGRR) DAREAR=0D0 + IF (RK.LE.0D0) THEN + AREA=0D0 + RETURN + ELSE +C +C +C PRESET ACCESSIBLE SOLID ANGLES TO "FREE" SPHERES. +C + AREA=FOURPI + ENDIF + 10 NCROSS=0 + EPSK=EPSI*RK +C +C CHECK IF SPHERE K NOT BURIED; SET UP FLAG FOR OVERLAP OF K. +C ------------------------------------------------------------- +C + DO 20 I=1,NAT + IF (I.EQ.K.OR.RAD(I).LE.0D0) GOTO 20 + IF (RK+RAD(I)-RLIO(IJ0(I,K)).LT.EPSK) GOTO 20 + IF (RLIO(IJ0(I,K))-ABS(RK-RAD(I)).LT.EPSK) THEN + IF (RK.LE.RAD(I)) THEN +C SPHERE K IMBEDDED IN SPHERE I + AREA=0D0 + RETURN + ENDIF + ELSE + NCROSS=NCROSS+1 + NC(NCROSS) = I + NCNCT(NCROSS,1)=I + ENDIF + 20 CONTINUE + IF (NCROSS.EQ.0) THEN +C THE SPHERE K IS FREE. + RETURN + ENDIF +C +C SET UP DATA FOR SPHERICAL SEGMENTS (SS). +C ---------------------------------------- +C + RKINV=0.5D0/RK + RK2=RK**2 + DO 23 I=1,NCROSS + LI=NCNCT(I,1) + RIKINV=1.D0/RLIO(IJ0(LI,K)) +C COSINE OF HALF-ANGLE OF CONE. + CTHETA(I,I)=RKINV*(RLIO(IJ0(LI,K))+(RK2-RAD(LI)**2)*RIKINV) + STHETA(I)=SQRT(1.D0-CTHETA(I,I)**2) + CONECT(I,I)=.FALSE. +C DIRECTOR COSINES OF VECTOR KI. + COSN(1,I,I)=URLIO(1,K,LI) + COSN(2,I,I)=URLIO(2,K,LI) + COSN(3,I,I)=URLIO(3,K,LI) + IF (LGRX.OR.LGRR) THEN + X=-CTHETA(I,I)/STHETA(I) + IF (LGRX) THEN + DRCTHT=RKINV*(1D0-(RK2-RAD(LI)**2)*RIKINV**2) + DO 22 ICOR=1,3 + DCTETA(ICOR,I,I)=DRCTHT*COSN(ICOR,I,I) + DSTETA(ICOR,I)=X*DCTETA(ICOR,I,I) + COSNI=COSN(ICOR,I,I)*RIKINV + DO 21 JCOR=1,ICOR + COSNIJ=-COSNI*COSN(JCOR,I,I) + DCOSN(ICOR,JCOR,I,I)=COSNIJ + 21 DCOSN(JCOR,ICOR,I,I)=COSNIJ + 22 DCOSN(ICOR,ICOR,I,I)=DCOSN(ICOR,ICOR,I,I)+RIKINV + ENDIF + IF (LGRR) THEN + DCTETR(I)=RIKINV-CTHETA(I,I)/RK + DSTETR(I)=X*DCTETR(I) + ENDIF + ENDIF + 23 CONTINUE +C +C CONNECTIVITY MATRIX BETWEEN SS. +C ------------------------------- +C TRUE IIF CONNECTED. DIAGONAL=TRUE IIF IMBEDDED. +C + DO 40 I=2,NCROSS + DO 30 J=1,I-1 + IF (CONECT(J,J)) GOTO 30 + CISJ=CTHETA(I,I)*STHETA(J) + SICJ=STHETA(I)*CTHETA(J,J) + SISJ=STHETA(I)*STHETA(J) +C COSINES OF ANGLES BETWEEN VECTOR KI AND KJ. + CTHETA(J,I)=DOTMN(COSN(1,I,I),COSN(1,J,J),3) + IF (LGRX) THEN + DO 24 ICOR=1,3 + DCTETA(ICOR,I,J)=DOTMN(DCOSN(1,ICOR,I,I),COSN(1,J,J),3) + 24 DCTETA(ICOR,J,I)=DOTMN(COSN(1,I,I),DCOSN(1,ICOR,J,J),3) + ENDIF + TIJ=CTHETA(J,I)-CTHETA(I,I)*CTHETA(J,J) + IF (TIJ.GT.SISJ-EPSI*ABS(CISJ-SICJ)) THEN + IF (CTHETA(J,J).GT.CTHETA(I,I)) THEN +C SS J IMBEDDED IN SS I + CONECT(J,J)=.TRUE. + ELSE +C SS I IMBEDDED IN SS J + CONECT(I,I)=.TRUE. + GOTO 40 + ENDIF + ELSE + EPSIJ=EPSI*(SICJ+CISJ) + IF (SICJ+CISJ.GE.0D0) THEN + CONECT(J,I)=TIJ.GT.EPSIJ-SISJ + ELSE + IF (TIJ.LE.-SISJ-EPSIJ) THEN +C SPHERE K FULLY COVERED BY SS I PLUS SS J + NCROSS=0 + AREA=0D0 + RETURN + ELSE + CONECT(J,I)=.TRUE. + ENDIF + ENDIF + CONECT(I,J)=CONECT(J,I) + ENDIF + 30 CONTINUE + 40 CONTINUE +C +C SUM THE CONTRIBUTIONS FROM ISOLATED SS. +C --------------------------------------- +C + ASLICE=0D0 + IF (LGRX) CALL SCOPYMN(3*(NCROSS+1),0D0,0,DCASLC(1,0),1) + IF (LGRR) DASLIR=0D0 + NCLUST=0 + DO 60 I=1,NCROSS + IF (CONECT(I,I)) GOTO 60 + DO 50 J=1,NCROSS + IF (CONECT(J,J)) GOTO 50 + IF (CONECT(J,I)) THEN + NCLUST=NCLUST+1 +C LABEL OF NON-ISOLATED SS: LAB + LAB(NCLUST)=I +C PREPARE DATA FOR FURTHER USE + WORK(NCLUST)=CTHETA(I,I)+EPSI*STHETA(I) + WORK(NCLUST+NCROSS)=CTHETA(I,I)-EPSI*STHETA(I) + GOTO 60 + ENDIF + 50 CONTINUE + ASLICE=ASLICE+1.0D0-CTHETA(I,I) + IF (LGRX) THEN + DO 55 ICOR=1,3 + DCASLC(ICOR,I)=DCASLC(ICOR,I)-DCTETA(ICOR,I,I) + 55 DCASLC(ICOR,0)=DCASLC(ICOR,0)+DCTETA(ICOR,I,I) + ENDIF + IF (LGRR) DASLIR=DASLIR-DCTETR(I) + 60 CONTINUE + ASLICE=TWOPI*ASLICE + IF (LGRX) CALL DSCALMN (3*(NCROSS+1),TWOPI,DCASLC(1,0),1) + IF (LGRR) DASLIR=DASLIR*TWOPI + IF (NCLUST.EQ.0) THEN +C NO CLUSTER OF SS REMAINS ON SPHERE K. + AREA=FOURPI-ASLICE + IF (LGRX) THEN + DO 61 I=0,NCROSS + DO 61 ICOR=1,3 + 61 DAREA(ICOR,I)=-DCASLC(ICOR,I) + ENDIF + IF (LGRR) DAREAR=-DASLIR + RETURN + ENDIF +C +C "FREE" INTERSECTIONS BETWEEN CLUSTERED SS. +C ------------------------------------------ +C + NFREE=0 + NCNCT(MXSS,LAB(1))=0 + DO 80 I=2,NCLUST + LI=LAB(I) + NCNCT(MXSS,LI)=0 + DO 80 J=1,I-1 + LJ=LAB(J) + IF (CONECT(LJ,LI)) THEN +C DIRECTOR COSINES OF INTERSECTION POINTS BETWEEN SS + SIN2IJ=1.D0/(1.D0-CTHETA(LJ,LI)**2) + AIJ=(CTHETA(LI,LI)-CTHETA(LJ,LJ)*CTHETA(LJ,LI))*SIN2IJ + BIJ=(CTHETA(LJ,LJ)-CTHETA(LI,LI)*CTHETA(LJ,LI))*SIN2IJ + CIJ=SQRT((1.D0-AIJ*CTHETA(LI,LI)-BIJ*CTHETA(LJ,LJ))*SIN2IJ) +C COMPUTE VN = COSN(LI,LI) X COSN(LJ,LJ). + VN(1)=COSN(2,LI,LI)*COSN(3,LJ,LJ)-COSN(3,LI,LI)*COSN(2,LJ,LJ) + VN(2)=COSN(3,LI,LI)*COSN(1,LJ,LJ)-COSN(1,LI,LI)*COSN(3,LJ,LJ) + VN(3)=COSN(1,LI,LI)*COSN(2,LJ,LJ)-COSN(2,LI,LI)*COSN(1,LJ,LJ) + IF (LGRX.OR.LGRR) THEN + DCIJ=0.5D0/CIJ + IF (LGRX) THEN + DO 62 ICOR=1,3 + DIVN(1,ICOR)=DCOSN(2,ICOR,LI,LI)*COSN(3,LJ,LJ) + 1 -DCOSN(3,ICOR,LI,LI)*COSN(2,LJ,LJ) + DIVN(2,ICOR)=DCOSN(3,ICOR,LI,LI)*COSN(1,LJ,LJ) + 1 -DCOSN(1,ICOR,LI,LI)*COSN(3,LJ,LJ) + DIVN(3,ICOR)=DCOSN(1,ICOR,LI,LI)*COSN(2,LJ,LJ) + 1 -DCOSN(2,ICOR,LI,LI)*COSN(1,LJ,LJ) + DJVN(1,ICOR)=COSN(2,LI,LI)*DCOSN(3,ICOR,LJ,LJ) + 1 -COSN(3,LI,LI)*DCOSN(2,ICOR,LJ,LJ) + DJVN(2,ICOR)=COSN(3,LI,LI)*DCOSN(1,ICOR,LJ,LJ) + 1 -COSN(1,LI,LI)*DCOSN(3,ICOR,LJ,LJ) + 62 DJVN(3,ICOR)=COSN(1,LI,LI)*DCOSN(2,ICOR,LJ,LJ) + 1 -COSN(2,LI,LI)*DCOSN(1,ICOR,LJ,LJ) + DSN2IJ=2.0D0*CTHETA(LJ,LI)*SIN2IJ + ENDIF + IF (LGRR) THEN + DAIJR=(DCTETR(LI)-DCTETR(LJ)*CTHETA(LJ,LI))*SIN2IJ + DBIJR=(DCTETR(LJ)-DCTETR(LI)*CTHETA(LJ,LI))*SIN2IJ + DCIJR=-SIN2IJ*(DAIJR*CTHETA(LI,LI)+AIJ*DCTETR(LI) + 1 +DBIJR*CTHETA(LJ,LJ)+BIJ*DCTETR(LJ))*DCIJ + ENDIF + ENDIF +C + DO 63 ICOR=1,3 +C +C ----- INTERSECTION PIJ ----- +C + COSN(ICOR,LI,LJ)=AIJ*COSN(ICOR,LI,LI)+BIJ*COSN(ICOR,LJ,LJ) + 1 +CIJ*VN(ICOR) +C +C ----- INTERSECTION PJI ----- +C + COSN(ICOR,LJ,LI)=AIJ*COSN(ICOR,LI,LI)+BIJ*COSN(ICOR,LJ,LJ) + 1 -CIJ*VN(ICOR) + IF (LGRX) THEN + DIS2IJ=DSN2IJ*DCTETA(ICOR,LI,LJ) + DJS2IJ=DSN2IJ*DCTETA(ICOR,LJ,LI) + DIAIJ(ICOR) = + 1 (DCTETA(ICOR,LI,LI)-CTHETA(LJ,LJ)*DCTETA(ICOR,LI,LJ)) + 2 *SIN2IJ+AIJ*DIS2IJ + DJAIJ(ICOR) = + 1 (-DCTETA(ICOR,LJ,LJ)*CTHETA(LJ,LI)-CTHETA(LJ,LJ) + 2 *DCTETA(ICOR,LJ,LI))*SIN2IJ + AIJ*DJS2IJ + DJBIJ(ICOR) = + 1 (DCTETA(ICOR,LJ,LJ)-CTHETA(LI,LI)*DCTETA(ICOR,LJ,LI)) + 2 *SIN2IJ + BIJ*DJS2IJ + DIBIJ(ICOR) = + 1 (-DCTETA(ICOR,LI,LI)*CTHETA(LJ,LI)-CTHETA(LI,LI) + 2 *DCTETA(ICOR,LI,LJ))*SIN2IJ + BIJ*DIS2IJ + DICIJ(ICOR) = + 1 -DCIJ*((DIAIJ(ICOR)*CTHETA(LI,LI)+AIJ*DCTETA(ICOR,LI,LI) + 2 +DIBIJ(ICOR)*CTHETA(LJ,LJ))*SIN2IJ)+0.5D0*CIJ*DIS2IJ + DJCIJ(ICOR) = + 1 -DCIJ*((DJBIJ(ICOR)*CTHETA(LJ,LJ)+BIJ*DCTETA(ICOR,LJ,LJ) + 2 +DJAIJ(ICOR)*CTHETA(LI,LI))*SIN2IJ)+0.5D0*CIJ*DJS2IJ + ENDIF + IF (LGRR) THEN + DCOSNR(ICOR,LI,LJ)=DAIJR*COSN(ICOR,LI,LI) + 1 +DBIJR*COSN(ICOR,LJ,LJ)+DCIJR*VN(ICOR) + DCOSNR(ICOR,LJ,LI)=DAIJR*COSN(ICOR,LI,LI) + 1 +DBIJR*COSN(ICOR,LJ,LJ)-DCIJR*VN(ICOR) + ENDIF + 63 CONTINUE +C + IF (LGRX) THEN + DO 64 ICOR=1,3 + DO 64 JCOR=1,3 + DICOS(ICOR,JCOR)=DIAIJ(JCOR)*COSN(ICOR,LI,LI) + + 1 AIJ*DCOSN(ICOR,JCOR,LI,LI)+DIBIJ(JCOR)*COSN(ICOR,LJ,LJ) + DJCOS(ICOR,JCOR)=DJAIJ(JCOR)*COSN(ICOR,LI,LI) + + 1 DJBIJ(JCOR)*COSN(ICOR,LJ,LJ)+BIJ*DCOSN(ICOR,JCOR,LJ,LJ) + DIWIJ(ICOR,JCOR)=DICIJ(JCOR)*VN(ICOR)+CIJ*DIVN(ICOR,JCOR) + DJWIJ(ICOR,JCOR)=DJCIJ(JCOR)*VN(ICOR)+CIJ*DJVN(ICOR,JCOR) + DICOSN(ICOR,JCOR,LI,LJ)=DICOS(ICOR,JCOR)+DIWIJ(ICOR,JCOR) + DJCOSN(ICOR,JCOR,LI,LJ)=DJCOS(ICOR,JCOR)+DJWIJ(ICOR,JCOR) + DICOSN(ICOR,JCOR,LJ,LI)=DICOS(ICOR,JCOR)-DIWIJ(ICOR,JCOR) + 64 DJCOSN(ICOR,JCOR,LJ,LI)=DJCOS(ICOR,JCOR)-DJWIJ(ICOR,JCOR) + ENDIF +C +C ----- INTERSECTION PIJ ------ +C + FREEIJ=.TRUE. +C +C ----- INTERSECTION PJI ----- +C + FREEJI=.TRUE. +C +C ----- SET UP TABLE OF FREE INTERSECTIONS BETWEEN SS. ----- +C ----- HERE IS THE L**3 PART OF THE CODE ----- +C ----- (L= MEAN VALUE OF NEIGHBORS) ----- +C + DO 70 L=1,NCLUST + IF (L.EQ.I.OR.L.EQ.J) GOTO 70 + LL=LAB(L) + IF (CONECT(LL,LI).AND.CONECT(LL,LJ)) THEN + IF (FREEJI) THEN + CHEK=DOTMN(COSN(1,LJ,LI),COSN(1,LL,LL),3) + IF (CHEK.GT.WORK(L)) THEN + FREEJI=.FALSE. + ELSE + PCMIN=MIN(PCMIN,ABS(WORK(L+NCROSS)-CHEK)) + IF (CHEK.GE.WORK(L+NCROSS)) THEN +C +C ----- SPHERES K, LAB(I), LAB(J), LAB(L) SHARE A POINT ----- +C ----- (AT THRESHOLD EPSI IN ANGLE). INCREASE THE ----- +C ----- RADIUS OF SPHERE K AND RESTART STUDY OF SPHERE K ----- +C + RK=RK*(1D0+4D0*EPSI) + JPCNT=JPCNT + 1 + GOTO 10 + ENDIF + ENDIF + ENDIF + IF (FREEIJ) THEN + CHEK=DOTMN(COSN(1,LI,LJ),COSN(1,LL,LL),3) + IF (CHEK.GT.WORK(L)) THEN + FREEIJ=.FALSE. + ELSE IF (CHEK.GE.WORK(L+NCROSS)) THEN + RK=RK*(1D0+4D0*EPSI) + GOTO 10 + ENDIF + ENDIF + IF (.NOT.(FREEIJ.OR.FREEJI)) GOTO 80 + ENDIF + 70 CONTINUE +C +C ----- ONE OR TWO FREE INTERSECTIONS FOUND. UPDATE TABLE ----- +C + CTHETA(LI,LJ)=CTHETA(LJ,LI) + IF (FREEJI) THEN + NFREE=NFREE+1 + M=NCNCT(MXSS,LI)+1 + NCNCT(M,LI)=LJ + NCNCT(MXSS,LI)=M + M=NCNCT(MXSS,LJ)+1 + NCNCT(M,LJ)=-LI + NCNCT(MXSS,LJ)=M + ENDIF + IF (FREEIJ) THEN + NFREE=NFREE+1 + M=NCNCT(MXSS,LI)+1 + NCNCT(M,LI)=-LJ + NCNCT(MXSS,LI)=M + M=NCNCT(MXSS,LJ)+1 + NCNCT(M,LJ)=LI + NCNCT(MXSS,LJ)=M + ENDIF + ENDIF + 80 CONTINUE + IF (NFREE.EQ.0) THEN +C SPHERE K BURIED BY A CLUSTER OF SS + NCROSS = 0 + AREA=0D0 + RETURN + ENDIF +C +C ORIENTED ANGLES OF SLICES OF SS. +C -------------------------------- +C + APOLY=0D0 + IF (LGRX) CALL SCOPYMN (3*(NCROSS+1),0D0,0,DCAPLY(1,0),1) + IF (LGRR) DAPOLR=0D0 + DO 130 I=1,NCLUST + IF (LGRX) CALL SCOPYMN (3*(NCROSS+1),0D0,0,DCAODD(1,0),1) + LI=LAB(I) + NPHI=NCNCT(MXSS,LI) + IF (NPHI.EQ.0) GOTO 130 + LJ=NCNCT(1,LI) + IF (LJ.GT.0) THEN + DO 83 ICOR=1,3 + CNIJ(ICOR) = COSN(ICOR,LJ,LI) + IF (LGRX.AND.LI.GE.LJ) THEN + DO 81 JCOR=1,3 + DICNIJ(ICOR,JCOR)=DICOSN(ICOR,JCOR,LJ,LI) + 81 DJCNIJ(ICOR,JCOR)=DJCOSN(ICOR,JCOR,LJ,LI) + ELSE IF (LGRX) THEN + DO 82 JCOR=1,3 + DICNIJ(ICOR,JCOR)=DJCOSN(ICOR,JCOR,LJ,LI) + 82 DJCNIJ(ICOR,JCOR)=DICOSN(ICOR,JCOR,LJ,LI) + ENDIF + IF (LGRR) DCNIJR(ICOR)=DCOSNR(ICOR,LJ,LI) + 83 CONTINUE + ELSE + LJ=-LJ + DO 86 ICOR=1,3 + CNIJ(ICOR)=COSN(ICOR,LI,LJ) + IF (LGRX.AND.LI.GE.LJ) THEN + DO 84 JCOR=1,3 + DICNIJ(ICOR,JCOR)=DICOSN(ICOR,JCOR,LI,LJ) + 84 DJCNIJ(ICOR,JCOR)=DJCOSN(ICOR,JCOR,LI,LJ) + ELSE IF (LGRX) THEN + DO 85 JCOR=1,3 + DICNIJ(ICOR,JCOR)=DJCOSN(ICOR,JCOR,LI,LJ) + 85 DJCNIJ(ICOR,JCOR)=DICOSN(ICOR,JCOR,LI,LJ) + ENDIF + IF (LGRR) DCNIJR(ICOR)=DCOSNR(ICOR,LI,LJ) + 86 CONTINUE + NCNCT(1,LI)=LJ + ENDIF + SIT(LI)=1.0D0/STHETA(LI) + C2I=CTHETA(LI,LI)**2 + VIJ(1)=COSN(2,LI,LI)*CNIJ(3)-COSN(3,LI,LI)*CNIJ(2) + VIJ(2)=COSN(3,LI,LI)*CNIJ(1)-COSN(1,LI,LI)*CNIJ(3) + VIJ(3)=COSN(1,LI,LI)*CNIJ(2)-COSN(2,LI,LI)*CNIJ(1) + LPOLY =DOTMN(VIJ,COSN(1,LJ,LJ),3).GT.0D0 +C + IF (LGRX.OR.LGRR) THEN + X=2D0*CTHETA(LI,LI) + Y=-SIT(LI)**2 + IF (LGRX) THEN + DO 87 ICOR=1,3 + DCSIT(ICOR,LI)=Y*DSTETA(ICOR,LI) + DCC2I(ICOR)=X*DCTETA(ICOR,LI,LI) + DIVIJ(1,ICOR)= + 1 DCOSN(2,ICOR,LI,LI)*CNIJ(3)-DCOSN(3,ICOR,LI,LI)*CNIJ(2) + 2 +COSN(2,LI,LI)*DICNIJ(3,ICOR)-COSN(3,LI,LI)*DICNIJ(2,ICOR) + DIVIJ(2,ICOR)= + 1 DCOSN(3,ICOR,LI,LI)*CNIJ(1)-DCOSN(1,ICOR,LI,LI)*CNIJ(3) + 2 +COSN(3,LI,LI)*DICNIJ(1,ICOR)-COSN(1,LI,LI)*DICNIJ(3,ICOR) + DIVIJ(3,ICOR)= + 1 DCOSN(1,ICOR,LI,LI)*CNIJ(2)-DCOSN(2,ICOR,LI,LI)*CNIJ(1) + 2 +COSN(1,LI,LI)*DICNIJ(2,ICOR)-COSN(2,LI,LI)*DICNIJ(1,ICOR) + DJVIJ(1,ICOR)= + 1 COSN(2,LI,LI)*DJCNIJ(3,ICOR)-COSN(3,LI,LI)*DJCNIJ(2,ICOR) + DJVIJ(2,ICOR)= + 1 COSN(3,LI,LI)*DJCNIJ(1,ICOR)-COSN(1,LI,LI)*DJCNIJ(3,ICOR) + 87 DJVIJ(3,ICOR)= + 1 COSN(1,LI,LI)*DJCNIJ(2,ICOR)-COSN(2,LI,LI)*DJCNIJ(1,ICOR) + ENDIF + IF (LGRR) THEN + DSITR(LI)=Y*DSTETR(LI) + DC2IR=X*DCTETR(LI) + DVIJR(1)=COSN(2,LI,LI)*DCNIJR(3)-COSN(3,LI,LI)*DCNIJR(2) + DVIJR(2)=COSN(3,LI,LI)*DCNIJR(1)-COSN(1,LI,LI)*DCNIJR(3) + DVIJR(3)=COSN(1,LI,LI)*DCNIJR(2)-COSN(2,LI,LI)*DCNIJR(1) + ENDIF + ENDIF +C +C MAKE A LOOP OVER THE SS I, COMPUTE DIHEDRALS BETWEEN 0 AND TWOPI +C + DO 95 J=2,NPHI + LJ=NCNCT(J,LI) + IF (LJ.GT.0) THEN + DO 90 ICOR=1,3 + CNIK(ICOR) = COSN(ICOR,LJ,LI) + IF (LGRX.AND.LI.GE.LJ) THEN + DO 88 JCOR=1,3 + DICNIK(ICOR,JCOR)=DICOSN(ICOR,JCOR,LJ,LI) + 88 DKCNIK(ICOR,JCOR)=DJCOSN(ICOR,JCOR,LJ,LI) + ELSE IF (LGRX) THEN + DO 89 JCOR=1,3 + DICNIK(ICOR,JCOR)=DJCOSN(ICOR,JCOR,LJ,LI) + 89 DKCNIK(ICOR,JCOR)=DICOSN(ICOR,JCOR,LJ,LI) + ENDIF + IF (LGRR) DCNIKR(ICOR)=DCOSNR(ICOR,LJ,LI) + 90 CONTINUE + ELSE + LJ =-LJ + DO 93 ICOR=1,3 + CNIK(ICOR)=COSN(ICOR,LI,LJ) + IF (LGRX.AND.LI.GE.LJ) THEN + DO 91 JCOR=1,3 + DICNIK(ICOR,JCOR)=DICOSN(ICOR,JCOR,LI,LJ) + 91 DKCNIK(ICOR,JCOR)=DJCOSN(ICOR,JCOR,LI,LJ) + ELSE IF (LGRX) THEN + DO 92 JCOR=1,3 + DICNIK(ICOR,JCOR)=DJCOSN(ICOR,JCOR,LI,LJ) + 92 DKCNIK(ICOR,JCOR)=DICOSN(ICOR,JCOR,LI,LJ) + ENDIF + IF (LGRR) DCNIKR(ICOR)=DCOSNR(ICOR,LI,LJ) + 93 CONTINUE + NCNCT(J,LI)=LJ + ENDIF +C +C NOTE: ACCURACY IS CRUCIAL HERE (USE ATAN2, NOT ACOS + ASIN) +C + X=DOTMN(VIJ,CNIK,3) + Y=DOTMN(CNIK,CNIJ,3)-C2I + WORK(J-1)=ATAN2(X,Y) + IF (WORK(J-1).LE.0D0) WORK(J-1)=WORK(J-1)+TWOPI + IF (LGRX.OR.LGRR) THEN + DX=-X/(X**2+Y**2) + DY= Y/(X**2+Y**2) + IF (LGRX) THEN +C CARTESIAN DERIVATIVES OF WORK + DO 94 ICOR=1,3 + DIX=DOTMN(DIVIJ(1,ICOR),CNIK,3)+ + $ DOTMN(VIJ,DICNIK(1,ICOR),3) + DIY=DOTMN(DICNIK(1,ICOR),CNIJ,3)+ + $ DOTMN(CNIK,DICNIJ(1,ICOR),3)-DCC2I(ICOR) + DJX=DOTMN(DJVIJ(1,ICOR),CNIK,3) + DJY=DOTMN(CNIK,DJCNIJ(1,ICOR),3) + DKX=DOTMN(VIJ,DKCNIK(1,ICOR),3) + DKY=DOTMN(DKCNIK(1,ICOR),CNIJ,3) + DIWORK(ICOR,J-1)= DY*DIX+DX*DIY + DJWORK(ICOR,J-1)= DY*DJX+DX*DJY + DKWORK(ICOR,J-1)= DY*DKX+DX*DKY + 94 D0WORK(ICOR,J-1)=-DIWORK(ICOR,J-1)-DJWORK(ICOR,J-1) + 1 -DKWORK(ICOR,J-1) + ENDIF + IF (LGRR) THEN +C DERIVATIVE OF WORK VS RAD(K) + DXR=DOTMN(DVIJR,CNIK,3)+DOTMN(VIJ,DCNIKR,3) + DYR=DOTMN(DCNIKR,CNIJ,3)+DOTMN(CNIK,DCNIJR,3)-DC2IR + DWORKR(J-1)=DY*DXR+DX*DYR + ENDIF + ENDIF + 95 CONTINUE + IF (NPHI.EQ.2) THEN +C TRIVIAL CASE, NO SORT NEEDED + AODD=WORK(1) + IF (LGRX) THEN + LJ=NCNCT(1,LI) + LK=NCNCT(2,LI) + DO 96 ICOR=1,3 + DCAODD(ICOR,LI)=DCAODD(ICOR,LI)+DIWORK(ICOR,1) + DCAODD(ICOR,LJ)=DCAODD(ICOR,LJ)+DJWORK(ICOR,1) + DCAODD(ICOR,LK)=DCAODD(ICOR,LK)+DKWORK(ICOR,1) + 96 DCAODD(ICOR,0 )=DCAODD(ICOR,0 )+D0WORK(ICOR,1) + ENDIF + IF (LGRR) DAODDR=DWORKR(1) + ELSE +C +C ----- SORT THE DIHEDRALS OF SS NO I IN ASCENDING ORDER ----- +C + DO 100 J=1,NPHI-2 + DO 100 L=J+1,NPHI-1 + IF (WORK(J).GT.WORK(L)) THEN + BUF=WORK(L) + WORK(L)=WORK(J) + WORK(J)=BUF + M=NCNCT(L+1,LI) + NCNCT(L+1,LI)=NCNCT(J+1,LI) + NCNCT(J+1,LI)=M + IF (LGRX) THEN + CALL DSWAPMN(3,DIWORK(1,L),1,DIWORK(1,J),1) + CALL DSWAPMN(3,DJWORK(1,L),1,DJWORK(1,J),1) + CALL DSWAPMN(3,DKWORK(1,L),1,DKWORK(1,J),1) + CALL DSWAPMN(3,D0WORK(1,L),1,D0WORK(1,J),1) + ENDIF + IF (LGRR) THEN + BUF=DWORKR(L) + DWORKR(L)=DWORKR(J) + DWORKR(J)=BUF + ENDIF + ENDIF + 100 CONTINUE +C +C COMPUTE THE SUM OF DIHEDRALS OF ORDERED SLICES (EVEN IN NUMBER) +C --------------------------------------------------------------- +C + AODD=WORK(1) + IF (LGRX) THEN + LJ=NCNCT(1,LI) + LK=NCNCT(2,LI) + DO 101 ICOR=1,3 + DCAODD(ICOR,LI)=DCAODD(ICOR,LI)+DIWORK(ICOR,1) + DCAODD(ICOR,LJ)=DCAODD(ICOR,LJ)+DJWORK(ICOR,1) + DCAODD(ICOR,LK)=DCAODD(ICOR,LK)+DKWORK(ICOR,1) + 101 DCAODD(ICOR,0 )=DCAODD(ICOR,0 )+D0WORK(ICOR,1) + ENDIF + IF (LGRR) DAODDR=DWORKR(1) + DO 110 J=3,NPHI-1,2 + AODD=AODD+WORK(J)-WORK(J-1) + IF (LGRX) THEN + LK=NCNCT(J,LI) + LL=NCNCT(J+1,LI) + DO 102 ICOR=1,3 + DCAODD(ICOR,LI)=DCAODD(ICOR,LI) + 1 +DIWORK(ICOR,J )-DIWORK(ICOR,J-1) + DCAODD(ICOR,LJ)=DCAODD(ICOR,LJ) + 1 +DJWORK(ICOR,J )-DJWORK(ICOR,J-1) + DCAODD(ICOR,LK)=DCAODD(ICOR,LK)-DKWORK(ICOR,J-1) + DCAODD(ICOR,LL)=DCAODD(ICOR,LL)+DKWORK(ICOR,J) + 102 DCAODD(ICOR,0 )=DCAODD(ICOR,0 ) + 1 +D0WORK(ICOR,J )-D0WORK(ICOR,J-1) + ENDIF + IF (LGRR) DAODDR=DAODDR+DWORKR(J)-DWORKR(J-1) + 110 CONTINUE + ENDIF + AEVEN=TWOPI-AODD + X=1D0-CTHETA(LI,LI) + IF (LPOLY) THEN +C ODD DIHEDRALS ARE VERTICE OF POLYGONS + APOLY=APOLY+AODD + ASLICE=ASLICE+AEVEN*X + IF (LGRX) THEN + DO 111 LT=0,NCROSS + DCAPLY(1,LT)=DCAPLY(1,LT)+DCAODD(1,LT) + DCASLC(1,LT)=DCASLC(1,LT)-DCAODD(1,LT)*X + DCAPLY(2,LT)=DCAPLY(2,LT)+DCAODD(2,LT) + DCASLC(2,LT)=DCASLC(2,LT)-DCAODD(2,LT)*X + DCAPLY(3,LT)=DCAPLY(3,LT)+DCAODD(3,LT) + 111 DCASLC(3,LT)=DCASLC(3,LT)-DCAODD(3,LT)*X + DCASLC(1,LI)=DCASLC(1,LI)-AEVEN*DCTETA(1,LI,LI) + DCASLC(1,0 )=DCASLC(1,0 )+AEVEN*DCTETA(1,LI,LI) + DCASLC(2,LI)=DCASLC(2,LI)-AEVEN*DCTETA(2,LI,LI) + DCASLC(2,0 )=DCASLC(2,0 )+AEVEN*DCTETA(2,LI,LI) + DCASLC(3,LI)=DCASLC(3,LI)-AEVEN*DCTETA(3,LI,LI) + DCASLC(3,0 )=DCASLC(3,0 )+AEVEN*DCTETA(3,LI,LI) + ENDIF + IF (LGRR) THEN + DAPOLR=DAPOLR+DAODDR + DASLIR=DASLIR-DAODDR*X-AEVEN*DCTETR(LI) + ENDIF + ELSE +C EVEN DIHEDRALS ARE VERTICE OF POLYGONS + APOLY=APOLY+AEVEN + ASLICE=ASLICE+AODD*X + IF (LGRX) THEN + DO 112 LT=0,NCROSS + DCAPLY(1,LT)=DCAPLY(1,LT)-DCAODD(1,LT) + DCASLC(1,LT)=DCASLC(1,LT)+DCAODD(1,LT)*X + DCAPLY(2,LT)=DCAPLY(2,LT)-DCAODD(2,LT) + DCASLC(2,LT)=DCASLC(2,LT)+DCAODD(2,LT)*X + DCAPLY(3,LT)=DCAPLY(3,LT)-DCAODD(3,LT) + 112 DCASLC(3,LT)=DCASLC(3,LT)+DCAODD(3,LT)*X + DCASLC(1,LI)=DCASLC(1,LI)-AODD*DCTETA(1,LI,LI) + DCASLC(1,0 )=DCASLC(1,0 )+AODD*DCTETA(1,LI,LI) + DCASLC(2,LI)=DCASLC(2,LI)-AODD*DCTETA(2,LI,LI) + DCASLC(2,0 )=DCASLC(2,0 )+AODD*DCTETA(2,LI,LI) + DCASLC(3,LI)=DCASLC(3,LI)-AODD*DCTETA(3,LI,LI) + DCASLC(3,0 )=DCASLC(3,0 )+AODD*DCTETA(3,LI,LI) + ENDIF + IF (LGRR) THEN + DAPOLR=DAPOLR-DAODDR + DASLIR=DASLIR+DAODDR*X-AODD*DCTETR(LI) + ENDIF +C REORDER LABELS SO THAT POLYGON'S VERTICE APPEAR AT ODD RANKS + M=NCNCT(1,LI) + DO 120 J=1,NPHI-1 + 120 NCNCT(J,LI)=NCNCT(J+1,LI) + NCNCT(NPHI,LI)=M + ENDIF + 130 CONTINUE +C +C NOW TO COUNT THE NUMBER OF POLYGONS AND TO SUM INTERSECTION ANGLES +C ------------------------------------------------------------------ +C + NPOLY=0 + DO 170 I=1,NCLUST + LI=LAB(I) + DO 160 J=2,NCNCT(MXSS,LI),2 + IF (NCNCT(J,LI).NE.0) THEN +C DEFINE THE "FIRST" VERTEX OF A POLYGON + IA=NCNCT(J-1,LI) + IB=LI + NCNCT(J,LI)=0 + PHI1=CTHETA(IA,IB)-CTHETA(IA,IA)*CTHETA(IB,IB) + PHI2=SIT(IA)*SIT(IB) + PHI=ACOS(PHI1*PHI2) + APOLY=APOLY+PHI + IF (LGRX.OR.LGRR) THEN + S1NPHI=1D0/SIN(PHI) + PHI1=S1NPHI*PHI1 + PHI2=S1NPHI*PHI2 + IF (LGRX) THEN + DO 131 ICOR=1,3 + DAPHI(ICOR)=(DCTETA(ICOR,IA,IB)-CTHETA(IB,IB) + 1 *DCTETA(ICOR,IA,IA))*PHI2+PHI1*DCSIT(ICOR,IA)*SIT(IB) + DBPHI(ICOR)=(DCTETA(ICOR,IB,IA)-CTHETA(IA,IA) + 1 *DCTETA(ICOR,IB,IB))*PHI2+PHI1*SIT(IA)*DCSIT(ICOR,IB) + DCAPLY(ICOR,IA)=DCAPLY(ICOR,IA)-DAPHI(ICOR) + DCAPLY(ICOR,IB)=DCAPLY(ICOR,IB)-DBPHI(ICOR) + 131 DCAPLY(ICOR,0 )=DCAPLY(ICOR,0 )+DAPHI(ICOR)+DBPHI(ICOR) + ENDIF + IF (LGRR) THEN + DPHIR= + 1 (DCTETR(IA)*CTHETA(IB,IB)+CTHETA(IA,IA)*DCTETR(IB))*PHI2 + 2 -PHI1*(DSITR(IA)*SIT(IB)+SIT(IA)*DSITR(IB)) + DAPOLR=DAPOLR+DPHIR + ENDIF + ENDIF +C +C FOLLOW AND ERASE VERTICES OF THE CURRENT POLYGON +C + 140 L=2 + 141 IF (L.GT.NCNCT(MXSS,IA)) GOTO 150 + IF (NCNCT(L,IA).EQ.IB) THEN + IBOLD=IB + NCNCT(L,IA)=0 + IB=IA + IA=NCNCT(L-1,IB) + IF (IA.NE.IBOLD) THEN + PHI1=CTHETA(IA,IB)-CTHETA(IA,IA)*CTHETA(IB,IB) + PHI2=SIT(IA)*SIT(IB) + PHI=ACOS(PHI1*PHI2) + IF (LGRX.OR.LGRR) THEN + S1NPHI=1D0/SIN(PHI) + PHI1=S1NPHI*PHI1 + PHI2=S1NPHI*PHI2 + IF (LGRX) THEN + DO 142 ICOR=1,3 + DAPHI(ICOR)=(DCTETA(ICOR,IA,IB) + 1 -CTHETA(IB,IB)*DCTETA(ICOR,IA,IA))*PHI2 + 2 +PHI1*DCSIT(ICOR,IA)*SIT(IB) + 142 DBPHI(ICOR)=(DCTETA(ICOR,IB,IA) + 1 -CTHETA(IA,IA)*DCTETA(ICOR,IB,IB))*PHI2 + 2 +PHI1*SIT(IA)*DCSIT(ICOR,IB) + ENDIF + IF (LGRR) THEN + DPHIR=(DCTETR(IA)*CTHETA(IB,IB)+ + 1 CTHETA(IA,IA)*DCTETR(IB))*PHI2 + 2 -PHI1*(DSITR(IA)*SIT(IB)+SIT(IA)*DSITR(IB)) + ENDIF + ENDIF + ELSE + IB=LI + IA=NCNCT(J-1,LI) + ENDIF + APOLY=APOLY+PHI + IF (LGRX) THEN + DO 143 ICOR=1,3 + DCAPLY(ICOR,IA)=DCAPLY(ICOR,IA)-DAPHI(ICOR) + DCAPLY(ICOR,IB)=DCAPLY(ICOR,IB)-DBPHI(ICOR) + 143 DCAPLY(ICOR,0 )=DCAPLY(ICOR,0 )+DAPHI(ICOR)+DBPHI(ICOR) + ENDIF + IF (LGRR) DAPOLR=DAPOLR+DPHIR + GOTO 140 + ENDIF + L=L+2 + GOTO 141 +C ORIENTED LOOP OVER CURRENT POLYGON COMPLETED + 150 NPOLY=NPOLY+1 + ENDIF + 160 CONTINUE + 170 CONTINUE +C +C ACCESSIBLE SOLID ANGLE FOR SPHERE K. +C ------------------------------------ +C POLYGONS IMBEDDED IN OTHERS TAKEN INTO ACCOUNT BY MODULO(4PI) +C + APOLY=APOLY+(NPOLY-NFREE)*TWOPI + AREA=FOURPI-ASLICE-MOD(APOLY,FOURPI) + IF (LGRX) THEN + DO 200 I=0,NCROSS + DAREA(1,I)=-DCASLC(1,I)-DCAPLY(1,I) + DAREA(2,I)=-DCASLC(2,I)-DCAPLY(2,I) + 200 DAREA(3,I)=-DCASLC(3,I)-DCAPLY(3,I) + ENDIF + IF (LGRR) DAREAR=-DASLIR-DAPOLR + RETURN + END +C +C ************** LIBRARIES FOR DAREAL SUBROUTINE ************** +C +C * SCOPYMN +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C COPIES A VECTOR, X, TO A VECTOR, Y. +C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. +C JACK DONGARRA, LINPACK, 3/11/78. +C +C LINPACK ROUTINE DCOPY CONVERTED TO SCOPY BY GCL 03/25/92 +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + SUBROUTINE SCOPYMN(N,DX,INCX,DY,INCY) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION DX(1),DY(1) +C + IF(N.LE.0)RETURN + IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 +C +C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS +C NOT EQUAL TO 1 +C + IX = 1 + IY = 1 + IF(INCX.LT.0)IX = (-N+1)*INCX + 1 + IF(INCY.LT.0)IY = (-N+1)*INCY + 1 + DO 10 I = 1,N + DY(IY) = DX(IX) + IX = IX + INCX + IY = IY + INCY + 10 CONTINUE + RETURN +C +C CODE FOR BOTH INCREMENTS EQUAL TO 1 +C +C +C CLEAN-UP LOOP +C + 20 M = MOD(N,7) + IF( M .EQ. 0 ) GO TO 40 + DO 30 I = 1,M + DY(I) = DX(I) + 30 CONTINUE + IF( N .LT. 7 ) RETURN + 40 MP1 = M + 1 + DO 50 I = MP1,N,7 + DY(I) = DX(I) + DY(I + 1) = DX(I + 1) + DY(I + 2) = DX(I + 2) + DY(I + 3) = DX(I + 3) + DY(I + 4) = DX(I + 4) + DY(I + 5) = DX(I + 5) + DY(I + 6) = DX(I + 6) + 50 CONTINUE + RETURN + END +C +C * DSCALMN + SUBROUTINE DSCALMN(N,DA,DX,INCX) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) +C +C SCALES A VECTOR BY A CONSTANT. +C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. +C JACK DONGARRA, LINPACK, 3/11/78. + DIMENSION DX(*) + IF(N.LE.0)RETURN + IF(INCX.EQ.1)GO TO 20 +C +C CODE FOR INCREMENT NOT EQUAL TO 1 +C + NINCX = N*INCX + DO 10 I = 1,NINCX,INCX + DX(I) = DA*DX(I) + 10 CONTINUE + RETURN +C +C CODE FOR INCREMENT EQUAL TO 1 +C +C +C CLEAN-UP LOOP +C + 20 M = MOD(N,5) + IF( M .EQ. 0 ) GO TO 40 + DO 30 I = 1,M + DX(I) = DA*DX(I) + 30 CONTINUE + IF( N .LT. 5 ) RETURN + 40 MP1 = M + 1 + DO 50 I = MP1,N,5 + DX(I) = DA*DX(I) + DX(I + 1) = DA*DX(I + 1) + DX(I + 2) = DA*DX(I + 2) + DX(I + 3) = DA*DX(I + 3) + DX(I + 4) = DA*DX(I + 4) + 50 CONTINUE + RETURN + END +C +C * DOTMN +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C DOT PRODUCT +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + FUNCTION DOTMN(X,Y,N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION X(*), Y(*) + DOTMN=DDOTMN(N,X,1,Y,1) + RETURN + END +C +C * DDOTMN + DOUBLE PRECISION FUNCTION DDOTMN(N,DX,INCX,DY,INCY) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION DX(*),DY(*) +C Purpose +C ======= +C forms the dot product of two vectors. +C uses unrolled loops for increments equal to one. +C jack dongarra, linpack, 3/11/78. +C modified 12/3/93, array(1) declarations changed to array(*) + INTRINSIC MOD + DDOTMN = 0.0d0 + DTEMP = 0.0d0 + IF (N.LE.0) RETURN + IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 + IX = 1 + IY = 1 + IF (INCX.LT.0) IX = (-N+1)*INCX + 1 + IF (INCY.LT.0) IY = (-N+1)*INCY + 1 + DO 10 I = 1,N + DTEMP = DTEMP + DX(IX)*DY(IY) + IX = IX + INCX + IY = IY + INCY + 10 CONTINUE + DDOTMN = DTEMP + RETURN + 20 M = MOD(N,5) + IF (M.EQ.0) GO TO 40 + DO 30 I = 1,M + DTEMP = DTEMP + DX(I)*DY(I) + 30 CONTINUE + IF (N.LT.5) GO TO 60 + 40 MP1 = M + 1 + DO 50 I = MP1,N,5 + DTEMP = DTEMP + DX(I)*DY(I) + DX(I+1)*DY(I+1) + + + DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4) + 50 CONTINUE + 60 DDOTMN = DTEMP + RETURN + END +C +C * DSWAPMN + SUBROUTINE DSWAPMN(N,DX,INCX,DY,INCY) + IMPLICIT DOUBLE PRECISION(A-H,O-Z) +c +c interchanges two vectors. +c uses unrolled loops for increments equal one. +c jack dongarra, linpack, 3/11/78. +c modified 12/3/93, array(1) declarations changed to array(*) +c + DIMENSION dx(*),dy(*) +c + if(n.le.0)return + if(incx.eq.1.and.incy.eq.1)go to 20 +c +c code for unequal increments or equal increments not equal +c to 1 +c + ix = 1 + iy = 1 + if(incx.lt.0)ix = (-n+1)*incx + 1 + if(incy.lt.0)iy = (-n+1)*incy + 1 + do 10 i = 1,n + dtemp = dx(ix) + dx(ix) = dy(iy) + dy(iy) = dtemp + ix = ix + incx + iy = iy + incy + 10 continue + return +c +c code for both increments equal to 1 +c +c +c clean-up loop +c + 20 m = mod(n,3) + if( m .eq. 0 ) go to 40 + do 30 i = 1,m + dtemp = dx(i) + dx(i) = dy(i) + dy(i) = dtemp + 30 continue + if( n .lt. 3 ) return + 40 mp1 = m + 1 + do 50 i = mp1,n,3 + dtemp = dx(i) + dx(i) = dy(i) + dy(i) = dtemp + dtemp = dx(i + 1) + dx(i + 1) = dy(i + 1) + dy(i + 1) = dtemp + dtemp = dx(i + 2) + dx(i + 2) = dy(i + 2) + dy(i + 2) = dtemp + 50 continue + return + end +c <-- MN solvation models +c $Id$ diff --git a/gpu4pyscf/lib/solvent/mnsol_interface.f90 b/gpu4pyscf/lib/solvent/mnsol_interface.f90 new file mode 100644 index 00000000..4a24b65a --- /dev/null +++ b/gpu4pyscf/lib/solvent/mnsol_interface.f90 @@ -0,0 +1,18 @@ +subroutine mnsol_interface(c,atnum,nat,sola,solb,solc,solg,solh,soln,icds,gcds,areacds,dcds) + IMPLICIT none + integer(kind=4), intent(in) :: nat,icds + real(kind=8), intent(in) :: sola,solb,solc,solg,solh,soln + real(kind=8), intent(in) :: c(3,nat) + real(kind=8), intent(out):: dcds(3,nat) + integer(kind=4), intent(in) :: atnum(nat) + real(kind=8), intent(out) :: gcds, areacds + + real(kind=8), allocatable :: x(:) + + integer :: ixmem + call mnsol_xmem(nat, ixmem) + allocate(x(ixmem)) + call cdsset(icds,gcds,areacds,nat,c,atnum,dcds,x,sola,solb,solc,solg,solh,soln) + deallocate(x) + +end subroutine mnsol_interface diff --git a/gpu4pyscf/lib/solvent/mnsol_mem.F b/gpu4pyscf/lib/solvent/mnsol_mem.F new file mode 100644 index 00000000..029ee800 --- /dev/null +++ b/gpu4pyscf/lib/solvent/mnsol_mem.F @@ -0,0 +1,98 @@ +c copied from https://github.com/nwchemgit/nwchem/blob/master/src/solvation/mnsol_mem.F +c +c +c memory for the x array used in mnsol code +c + subroutine mnsol_xmem(nat,ineed) +c + implicit none +c + integer nat +c + integer LRAD + integer LCDSA + integer LAREA + integer LDATAR + integer LSTS + integer LCOT + integer LDSTS + integer LDCOTDR + integer LNC + integer LDAREA + integer LRLIO + integer LURLIO + integer LLAB + integer LNCNCT + integer LCONECT + integer LCTHETA + integer LSTHETA + integer LSIT + integer LDCSIT + integer LDJCOSN + integer LCOSN + integer LDSTETA + integer LDCTETA + integer LDCOSN + integer LWORK + integer LDIWORK + integer LDJWORK + integer LDKWORK + integer LD0WORK + integer LDWORKR + integer LDCAODD + integer LDSITR + integer LDCAPLY + integer LDICOSN + integer LDCASLC + integer LDCTETR + integer LDSTETR + integer LDCOSNR + integer LEND + integer INEED +c +c calculate memory +c + LRAD = 1 + LCDSA = LRAD + NAT + LAREA = LCDSA + NAT + LDATAR = LAREA + NAT + LSTS = LDATAR + 3 * NAT * NAT + LCOT = LSTS + NAT + LDSTS = LCOT + NAT * (NAT + 1) / 2 + LDCOTDR = LDSTS + 3 * NAT * NAT + LNC = LDCOTDR + NAT * (NAT + 1) / 2 + LDAREA = LNC + NAT + 1 + LRLIO = LDAREA + 3 * (NAT + 1) + LURLIO = LRLIO + NAT * (NAT + 1) / 2 + LLAB = LURLIO + 3 * NAT * NAT + LNCNCT = LLAB + NAT + LCONECT = LNCNCT + NAT * (2 * NAT + 1) + LCTHETA = LCONECT + NAT * NAT + LSTHETA = LCTHETA + NAT * NAT + LSIT = LSTHETA + NAT + LDCSIT = LSIT + NAT + LDJCOSN = LDCSIT + 3 * NAT + LCOSN = LDJCOSN + 3 * 3 * NAT * NAT + LDSTETA = LCOSN + 3 * NAT * NAT + LDCTETA = LDSTETA + 3 * NAT + LDCOSN = LDCTETA + 3 * NAT * NAT + LWORK = LDCOSN + 3 * 3 * NAT * NAT + LDIWORK = LWORK + 2 * NAT + 1 + LDJWORK = LDIWORK + 3 * NAT + LDKWORK = LDJWORK + 3 * NAT + LD0WORK = LDKWORK + 3 * NAT + LDWORKR = LD0WORK + 3 * NAT + LDCAODD = LDWORKR + NAT + LDSITR = LDCAODD + 3 * (NAT + 1) + LDCAPLY = LDSITR + NAT + LDICOSN = LDCAPLY + 3 * (NAT + 1) + LDCASLC = LDICOSN + 3 * 3 * NAT * NAT + LDCTETR = LDCASLC + 3 * (NAT + 1) + LDSTETR = LDCTETR + NAT + LDCOSNR = LDSTETR + NAT + LEND = LDCOSNR + 3 * NAT * NAT + INEED = LEND - LRAD +c + return + end +c $Id$ diff --git a/gpu4pyscf/lib/tests/test_cupy_helper.py b/gpu4pyscf/lib/tests/test_cupy_helper.py index 4ee8d3cd..af54dc0d 100644 --- a/gpu4pyscf/lib/tests/test_cupy_helper.py +++ b/gpu4pyscf/lib/tests/test_cupy_helper.py @@ -176,6 +176,22 @@ def test_cart2sph(self): a_sph1 = cart2sph(a_cart, axis=1, ang=4) assert cupy.linalg.norm(a_sph0 - a_sph1) < 1e-8 + a_cart = cupy.random.rand(10,21,11) + a_sph0 = cart2sph_cutensor(a_cart, axis=1, ang=5) + a_sph1 = cart2sph(a_cart, axis=1, ang=5) + assert cupy.linalg.norm(a_sph0 - a_sph1) < 1e-8 + + a_cart = cupy.random.rand(10,28,11) + a_sph0 = cart2sph_cutensor(a_cart, axis=1, ang=6) + a_sph1 = cart2sph(a_cart, axis=1, ang=6) + assert cupy.linalg.norm(a_sph0 - a_sph1) < 1e-8 + + ''' + a_cart = cupy.random.rand(10,36,11) + a_sph0 = cart2sph_cutensor(a_cart, axis=1, ang=7) + a_sph1 = cart2sph(a_cart, axis=1, ang=7) + assert cupy.linalg.norm(a_sph0 - a_sph1) < 1e-8 + ''' if __name__ == "__main__": print("Full tests for cupy helper module") unittest.main() diff --git a/gpu4pyscf/lib/tests/test_to_gpu.py b/gpu4pyscf/lib/tests/test_to_gpu.py index f55421c2..80f60b62 100644 --- a/gpu4pyscf/lib/tests/test_to_gpu.py +++ b/gpu4pyscf/lib/tests/test_to_gpu.py @@ -127,9 +127,6 @@ def test_df_RKS(self): h = hobj.kernel() assert numpy.abs(lib.fp(h) - 2.187025544697092) < 1e-4 - - - if __name__ == "__main__": print("Full tests for to_gpu module") unittest.main() \ No newline at end of file diff --git a/gpu4pyscf/solvent/_attach_solvent.py b/gpu4pyscf/solvent/_attach_solvent.py index d0bfe330..8e27bb54 100644 --- a/gpu4pyscf/solvent/_attach_solvent.py +++ b/gpu4pyscf/solvent/_attach_solvent.py @@ -16,12 +16,9 @@ import cupy from pyscf import lib from pyscf.lib import logger -from pyscf.solvent._attach_solvent import _Solvation from gpu4pyscf.lib.cupy_helper import tag_array from gpu4pyscf import scf -# NOTE: copied from pyscf, different from the latest version - def _for_scf(mf, solvent_obj, dm=None): '''Add solvent model to SCF (HF and DFT) method. @@ -33,120 +30,131 @@ def _for_scf(mf, solvent_obj, dm=None): mf.with_solvent = solvent_obj return mf - oldMF = mf.__class__ - if dm is not None: solvent_obj.e, solvent_obj.v = solvent_obj.kernel(dm) solvent_obj.frozen = True - class SCFWithSolvent(_Solvation, oldMF): - def __init__(self, mf, solvent): - self.__dict__.update(mf.__dict__) - self.with_solvent = solvent - self._keys.update(['with_solvent']) - - def dump_flags(self, verbose=None): - oldMF.dump_flags(self, verbose) - self.with_solvent.check_sanity() - self.with_solvent.dump_flags(verbose) - return self - - def reset(self, mol=None): - self.with_solvent.reset(mol) - return oldMF.reset(self, mol) - - # Note v_solvent should not be added to get_hcore for scf methods. - # get_hcore is overloaded by many post-HF methods. Modifying - # SCF.get_hcore may lead error. - - def get_veff(self, mol=None, dm=None, *args, **kwargs): - vhf = oldMF.get_veff(self, mol, dm, *args, **kwargs) - with_solvent = self.with_solvent - if not with_solvent.frozen: - with_solvent.e, with_solvent.v = with_solvent.kernel(dm) - e_solvent, v_solvent = with_solvent.e, with_solvent.v - - # NOTE: v_solvent should not be added to vhf in this place. This is - # because vhf is used as the reference for direct_scf in the next - # iteration. If v_solvent is added here, it may break direct SCF. - return tag_array(vhf, e_solvent=e_solvent, v_solvent=v_solvent) - - def get_fock(self, h1e=None, s1e=None, vhf=None, dm=None, cycle=-1, - diis=None, diis_start_cycle=None, - level_shift_factor=None, damp_factor=None): - # DIIS was called inside oldMF.get_fock. v_solvent, as a function of - # dm, should be extrapolated as well. To enable it, v_solvent has to be - # added to the fock matrix before DIIS was called. - if getattr(vhf, 'v_solvent', None) is None: - vhf = self.get_veff(self.mol, dm) - return oldMF.get_fock(self, h1e, s1e, vhf+vhf.v_solvent, dm, cycle, diis, - diis_start_cycle, level_shift_factor, damp_factor) - - def energy_elec(self, dm=None, h1e=None, vhf=None): - if dm is None: - dm = self.make_rdm1() - if getattr(vhf, 'e_solvent', None) is None: - vhf = self.get_veff(self.mol, dm) - - e_tot, e_coul = oldMF.energy_elec(self, dm, h1e, vhf) - e_solvent = vhf.e_solvent - if isinstance(e_solvent, cupy.ndarray): - e_solvent = e_solvent.get()[()] - e_tot += e_solvent - self.scf_summary['e_solvent'] = vhf.e_solvent.real - - if self.with_solvent.method.upper() == 'SMD': - if self.with_solvent.e_cds is None: - e_cds = self.with_solvent.get_cds() - self.with_solvent.e_cds = e_cds - else: - e_cds = self.with_solvent.e_cds - if isinstance(e_cds, cupy.ndarray): - e_cds = e_cds.get()[()] - e_tot += e_cds - self.scf_summary['e_cds'] = e_cds - logger.info(self, f'CDS correction = {e_cds:.15f}') - logger.info(self, 'Solvent Energy = %.15g', vhf.e_solvent) - - return e_tot, e_coul - - def nuc_grad_method(self): - grad_method = oldMF.nuc_grad_method(self) - return self.with_solvent.nuc_grad_method(grad_method) - - Gradients = nuc_grad_method - - def Hessian(self): - hess_method = oldMF.Hessian(self) - return self.with_solvent.Hessian(hess_method) - - def gen_response(self, *args, **kwargs): - vind = oldMF.gen_response(self, *args, **kwargs) - is_uhf = isinstance(self, scf.uhf.UHF) - # singlet=None is orbital hessian or CPHF type response function - singlet = kwargs.get('singlet', True) - singlet = singlet or singlet is None - def vind_with_solvent(dm1): - v = vind(dm1) - if self.with_solvent.equilibrium_solvation: - if is_uhf: - v_solvent = self.with_solvent._B_dot_x(dm1) - v += v_solvent[0] + v_solvent[1] - elif singlet: - v += self.with_solvent._B_dot_x(dm1) - return v - return vind_with_solvent - - def stability(self, *args, **kwargs): - # When computing orbital hessian, the second order derivatives of - # solvent energy needs to be computed. It is enabled by - # the attribute equilibrium_solvation in gen_response method. - # If solvent was frozen, its contribution is treated as the - # external potential. The response of solvent does not need to - # be considered in stability analysis. - with lib.temporary_env(self.with_solvent, - equilibrium_solvation=not self.with_solvent.frozen): - return oldMF.stability(self, *args, **kwargs) - - mf1 = SCFWithSolvent(mf, solvent_obj) - return mf1 \ No newline at end of file + sol_mf = SCFWithSolvent(mf, solvent_obj) + name = solvent_obj.__class__.__name__ + mf.__class__.__name__ + return lib.set_class(sol_mf, (SCFWithSolvent, mf.__class__), name) + +# 1. A tag to label the derived method class +class _Solvation: + pass + +class SCFWithSolvent(_Solvation): + from gpu4pyscf.lib.utils import to_gpu, device + + _keys = {'with_solvent'} + + def __init__(self, mf, solvent): + self.__dict__.update(mf.__dict__) + self.with_solvent = solvent + + def undo_solvent(self): + cls = self.__class__ + name_mixin = self.with_solvent.__class__.__name__ + obj = lib.view(self, lib.drop_class(cls, SCFWithSolvent, name_mixin)) + del obj.with_solvent + return obj + + def to_cpu(self): + from pyscf.solvent import _attach_solvent + solvent_obj = self.with_solvent.to_cpu() + obj = _attach_solvent._for_scf(self.undo_solvent().to_cpu(), solvent_obj) + return obj + + def dump_flags(self, verbose=None): + super().dump_flags(verbose) + self.with_solvent.check_sanity() + self.with_solvent.dump_flags(verbose) + return self + + def reset(self, mol=None): + self.with_solvent.reset(mol) + return super().reset(mol) + + def get_veff(self, mol=None, dm=None, *args, **kwargs): + vhf = super().get_veff(mol, dm, *args, **kwargs) + with_solvent = self.with_solvent + if not with_solvent.frozen: + with_solvent.e, with_solvent.v = with_solvent.kernel(dm) + e_solvent, v_solvent = with_solvent.e, with_solvent.v + return tag_array(vhf, e_solvent=e_solvent, v_solvent=v_solvent) + + def get_fock(self, h1e=None, s1e=None, vhf=None, dm=None, cycle=-1, + diis=None, diis_start_cycle=None, + level_shift_factor=None, damp_factor=None, fock_last=None): + # DIIS was called inside oldMF.get_fock. v_solvent, as a function of + # dm, should be extrapolated as well. To enable it, v_solvent has to be + # added to the fock matrix before DIIS was called. + if getattr(vhf, 'v_solvent', None) is None: + vhf = self.get_veff(self.mol, dm) + return super().get_fock(h1e, s1e, vhf+vhf.v_solvent, dm, cycle, diis, + diis_start_cycle, level_shift_factor, damp_factor) + + def energy_elec(self, dm=None, h1e=None, vhf=None): + if dm is None: + dm = self.make_rdm1() + if getattr(vhf, 'e_solvent', None) is None: + vhf = self.get_veff(self.mol, dm) + + e_tot, e_coul = super().energy_elec(dm, h1e, vhf) + e_solvent = vhf.e_solvent + if isinstance(e_solvent, cupy.ndarray): + e_solvent = e_solvent.get()[()] + e_tot += e_solvent + self.scf_summary['e_solvent'] = vhf.e_solvent.real + + if (hasattr(self.with_solvent, 'method') and self.with_solvent.method.upper() == 'SMD'): + if self.with_solvent.e_cds is None: + e_cds = self.with_solvent.get_cds() + self.with_solvent.e_cds = e_cds + else: + e_cds = self.with_solvent.e_cds + if isinstance(e_cds, cupy.ndarray): + e_cds = e_cds.get()[()] + e_tot += e_cds + self.scf_summary['e_cds'] = e_cds + logger.info(self, f'CDS correction = {e_cds:.15f}') + logger.info(self, 'Solvent Energy = %.15g', vhf.e_solvent) + + return e_tot, e_coul + + def nuc_grad_method(self): + grad_method = super().nuc_grad_method() + return self.with_solvent.nuc_grad_method(grad_method) + + Gradients = nuc_grad_method + + def Hessian(self): + hess_method = super().Hessian() + return self.with_solvent.Hessian(hess_method) + + def gen_response(self, *args, **kwargs): + vind = super().gen_response(*args, **kwargs) + is_uhf = isinstance(self, scf.uhf.UHF) + # singlet=None is orbital hessian or CPHF type response function + singlet = kwargs.get('singlet', True) + singlet = singlet or singlet is None + def vind_with_solvent(dm1): + v = vind(dm1) + if self.with_solvent.equilibrium_solvation: + if is_uhf: + v_solvent = self.with_solvent._B_dot_x(dm1) + v += v_solvent[0] + v_solvent[1] + elif singlet: + v += self.with_solvent._B_dot_x(dm1) + return v + return vind_with_solvent + + def stability(self, *args, **kwargs): + # When computing orbital hessian, the second order derivatives of + # solvent energy needs to be computed. It is enabled by + # the attribute equilibrium_solvation in gen_response method. + # If solvent was frozen, its contribution is treated as the + # external potential. The response of solvent does not need to + # be considered in stability analysis. + with lib.temporary_env(self.with_solvent, + equilibrium_solvation=not self.with_solvent.frozen): + return super().stability(*args, **kwargs) diff --git a/gpu4pyscf/solvent/grad/pcm.py b/gpu4pyscf/solvent/grad/pcm.py index d0f2576b..508daa31 100644 --- a/gpu4pyscf/solvent/grad/pcm.py +++ b/gpu4pyscf/solvent/grad/pcm.py @@ -181,7 +181,12 @@ def grad_nuc(pcmobj, dm): mol = pcmobj.mol log = logger.new_logger(mol, mol.verbose) t1 = log.init_timer() - if not pcmobj._intermediates or 'q_sym' not in pcmobj._intermediates: + if not pcmobj._intermediates: + pcmobj.build() + dm_cache = pcmobj._intermediates.get('dm', None) + if dm_cache is not None and cupy.linalg.norm(dm_cache - dm) < 1e-10: + pass + else: pcmobj._get_vind(dm) mol = pcmobj.mol @@ -215,17 +220,23 @@ def grad_qv(pcmobj, dm): ''' contributions due to integrals ''' - if not pcmobj._intermediates or 'q_sym' not in pcmobj._intermediates: + if not pcmobj._intermediates: + pcmobj.build() + dm_cache = pcmobj._intermediates.get('dm', None) + if dm_cache is not None and cupy.linalg.norm(dm_cache - dm) < 1e-10: + pass + else: pcmobj._get_vind(dm) mol = pcmobj.mol log = logger.new_logger(mol, mol.verbose) t1 = log.init_timer() - gridslice = pcmobj.surface['gslice_by_atom'] - q_sym = pcmobj._intermediates['q_sym'] + gridslice = pcmobj.surface['gslice_by_atom'] + charge_exp = pcmobj.surface['charge_exp'] + grid_coords = pcmobj.surface['grid_coords'] + q_sym = pcmobj._intermediates['q_sym'] - intopt = pcmobj.intopt - intopt.clear() - # rebuild with aosym + auxmol = gto.fakemol_for_charges(grid_coords.get(), expnt=charge_exp.get()**2) + intopt = int3c2e.VHFOpt(mol, auxmol, 'int2e') intopt.build(1e-14, diag_block_with_triu=True, aosym=False) coeff = intopt.coeff dm_cart = coeff @ dm @ coeff.T @@ -252,7 +263,12 @@ def grad_solver(pcmobj, dm): mol = pcmobj.mol log = logger.new_logger(mol, mol.verbose) t1 = log.init_timer() - if not pcmobj._intermediates or 'q_sym' not in pcmobj._intermediates: + if not pcmobj._intermediates: + pcmobj.build() + dm_cache = pcmobj._intermediates.get('dm', None) + if dm_cache is not None and cupy.linalg.norm(dm_cache - dm) < 1e-10: + pass + else: pcmobj._get_vind(dm) gridslice = pcmobj.surface['gslice_by_atom'] @@ -360,6 +376,8 @@ def make_grad_object(grad_method): (WithSolventGrad, grad_method.__class__), name) class WithSolventGrad: + from gpu4pyscf.lib.utils import to_gpu, device + _keys = {'de_solvent', 'de_solute'} def __init__(self, grad_method): @@ -375,6 +393,11 @@ def undo_solvent(self): del obj.de_solute return obj + def to_cpu(self): + from pyscf.solvent.grad import pcm # type: ignore + grad_method = self.undo_solvent().to_cpu() + return pcm.make_grad_object(grad_method) + def kernel(self, *args, dm=None, atmlst=None, **kwargs): dm = kwargs.pop('dm', None) if dm is None: diff --git a/gpu4pyscf/solvent/grad/smd.py b/gpu4pyscf/solvent/grad/smd.py index b8cc34ba..fe315fd6 100644 --- a/gpu4pyscf/solvent/grad/smd.py +++ b/gpu4pyscf/solvent/grad/smd.py @@ -14,7 +14,7 @@ # along with this program. If not, see . ''' -Gradient of PCM family solvent model +Gradient of SMD solvent model ''' # pylint: disable=C0103 @@ -23,207 +23,111 @@ #from cupyx import scipy, jit from pyscf import lib from pyscf.grad import rhf as rhf_grad -from pyscf.data import radii - from gpu4pyscf.solvent import pcm, smd from gpu4pyscf.solvent.grad import pcm as pcm_grad -from gpu4pyscf.solvent.smd import ( - sigma_water, sigma_n, sigma_alpha, sigma_beta, r_zz, swtich_function, - hartree2kcal) -from gpu4pyscf.lib.cupy_helper import contract from gpu4pyscf.lib import logger +from gpu4pyscf.lib.cupy_helper import contract -def grad_swtich_function(R, r, dr): - if R < r + dr: - return -np.exp(dr/(R-dr-r)) * dr / (R-dr-r)**2 - else: - return 0.0 +def get_cds(smdobj): + return smd.get_cds_legacy(smdobj)[1] -def atomic_surface_tension(symbols, coords, n, alpha, beta, water=True): +def grad_solver(smdobj, dm): ''' - - list of atomic symbols - - atomic coordinates in Anstrong - - solvent descriptors: n, alpha, beta + dE = 0.5*v* d(K^-1 R) *v + q*dv + v^T* d(K^-1 R)v = v^T*K^-1(dR - dK K^-1R)v = v^T K^-1(dR - dK q) ''' - - def get_bond_tension(bond): - if water: - return sigma_water.get(bond, 0.0) - t = 0.0 - t += sigma_n.get(bond, 0.0) * n - t += sigma_alpha.get(bond, 0.0) * alpha - t += sigma_beta.get(bond, 0.0) * beta - return t - - def get_atom_tension(sym_i): - if water: - return sigma_water.get(sym_i, 0.0) - t = 0.0 - t += sigma_n.get(sym_i, 0.0) * n - t += sigma_alpha.get(sym_i, 0.0) * alpha - t += sigma_beta.get(sym_i, 0.0) * beta - return t - natm = coords.shape[0] - ri_rj = coords[:,None,:] - coords[None,:,:] - rij = np.sum(ri_rj**2, axis=2)**0.5 - np.fill_diagonal(rij, 1) - drij = ri_rj / np.expand_dims(rij, axis=-1) - tensions = [] - for i, sym_i in enumerate(symbols): - if sym_i not in ['H', 'C', 'N', 'O', 'F', 'Si', 'S', 'Cl', 'Br']: - tensions.append(np.zeros([natm,3])) - continue - - tension = np.zeros([natm,3]) - if sym_i in ['F', 'Si', 'S', 'Cl', 'Br']: - tensions.append(tension) - continue - - if sym_i == 'H': - dt_HC = np.zeros([natm,3]) - dt_HO = np.zeros([natm,3]) - for j, sym_j in enumerate(symbols): - if sym_j == 'C': - r, dr = r_zz.get(('H','C'), (0.0, 0.0)) - dt_drij = grad_swtich_function(rij[i,j], r, dr) * drij[i,j] - dt_HC[i] += dt_drij - dt_HC[j] -= dt_drij - if sym_j == 'O': - r, dr = r_zz.get(('H','O'), (0.0, 0.0)) - dt_drij = grad_swtich_function(rij[i,j], r, dr) * drij[i,j] - dt_HO[i] += dt_drij - dt_HO[j] -= dt_drij - sig_HC = get_bond_tension(('H','C')) - sig_HO = get_bond_tension(('H','O')) - tension += sig_HC * dt_HC + sig_HO * dt_HO - tensions.append(tension) - continue - - if sym_i == 'C': - dt_CC = np.zeros([natm,3]) - dt_CN = np.zeros([natm,3]) - t_CN = 0.0 - for j, sym_j in enumerate(symbols): - if sym_j == 'C' and i != j: - r, dr = r_zz.get(('C', 'C'), (0.0, 0.0)) - dt_drij = grad_swtich_function(rij[i,j], r, dr) * drij[i,j] - dt_CC[i] += dt_drij - dt_CC[j] -= dt_drij - if sym_j == 'N': - r, dr = r_zz.get(('C', 'N'), (0.0, 0.0)) - t_CN += swtich_function(rij[i,j], r, dr) - dt_drij = grad_swtich_function(rij[i,j], r, dr) * drij[i,j] - dt_CN[i] += dt_drij - dt_CN[j] -= dt_drij - sig_CC = get_bond_tension(('C','C')) - sig_CN = get_bond_tension(('C','N')) - tension += sig_CC * dt_CC + sig_CN * (2 * t_CN * dt_CN) - tensions.append(tension) - continue - - if sym_i == 'N': - t_NC = 0.0 - dt_NC = np.zeros([natm,3]) - dt_NC3 = np.zeros([natm,3]) - for j, sym_j in enumerate(symbols): - if sym_j == 'C': - r, dr = r_zz.get(('N','C'), (0.0, 0.0)) - tk = 0.0 - dtk = np.zeros([natm,3]) - for k, sym_k in enumerate(symbols): - if k != i and k != j: - rjk, drjk = r_zz.get(('C', sym_k), (0.0, 0.0)) - tk += swtich_function(rij[j,k], rjk, drjk) - dtk_rjk = grad_swtich_function(rij[j,k], rjk, drjk) * drij[j,k] - dtk[j] += dtk_rjk - dtk[k] -= dtk_rjk - - dt_drij = grad_swtich_function(rij[i,j], r, dr) * drij[i,j] * tk**2 - dt_NC[i] += dt_drij - dt_NC[j] -= dt_drij - - t = swtich_function(rij[i,j], r, dr) - dt_NC += t * (2 * tk * dtk) - t_NC += t * tk**2 - - r, dr = r_zz.get(('N','C3'), (0.0, 0.0)) - dt_drij = grad_swtich_function(rij[i,j], r, dr) * drij[i,j] - dt_NC3[i] += dt_drij - dt_NC3[j] -= dt_drij - sig_NC = get_bond_tension(('N','C')) - sig_NC3= get_bond_tension(('N','C3')) - tension += sig_NC * (1.3 * t_NC**0.3 * dt_NC) + sig_NC3 * dt_NC3 - tensions.append(tension) - continue - - if sym_i == 'O': - dt_OC = np.zeros([natm,3]) - dt_ON = np.zeros([natm,3]) - dt_OO = np.zeros([natm,3]) - dt_OP = np.zeros([natm,3]) - for j, sym_j in enumerate(symbols): - if sym_j == 'C': - r, dr = r_zz.get(('O','C'), (0.0, 0.0)) - dt_drij = grad_swtich_function(rij[i,j], r, dr) * drij[i,j] - dt_OC[i] += dt_drij - dt_OC[j] -= dt_drij - if sym_j == 'N': - r, dr = r_zz.get(('O','N'), (0.0, 0.0)) - dt_drij = grad_swtich_function(rij[i,j], r, dr) * drij[i,j] - dt_ON[i] += dt_drij - dt_ON[j] -= dt_drij - if sym_j == 'O' and j != i: - r, dr = r_zz.get(('O','O'), (0.0, 0.0)) - dt_drij = grad_swtich_function(rij[i,j], r, dr) * drij[i,j] - dt_OO[i] += dt_drij - dt_OO[j] -= dt_drij - if sym_j == 'P': - r, dr = r_zz.get(('O','P'), (0.0, 0.0)) - dt_drij = grad_swtich_function(rij[i,j], r, dr) * drij[i,j] - dt_OP[i] += dt_drij - dt_OP[j] -= dt_drij - sig_OC = get_bond_tension(('O','C')) - sig_ON = get_bond_tension(('O','N')) - sig_OO = get_bond_tension(('O','O')) - sig_OP = get_bond_tension(('O','P')) - tension += sig_OC * dt_OC + sig_ON * dt_ON + sig_OO * dt_OO + sig_OP * dt_OP - tensions.append(tension) - continue - - return cupy.asarray(tensions) - -def get_cds(smdobj): mol = smdobj.mol - n, _, alpha, beta, gamma, _, phi, psi = smdobj.solvent_descriptors - symbols = [mol.atom_symbol(ia) for ia in range(mol.natm)] - coords = mol.atom_coords(unit='A') - if smdobj._solvent.lower() != 'water': - grad_tension = atomic_surface_tension(symbols, coords, n, alpha, beta, water=False) - atm_tension = smd.atomic_surface_tension(symbols, coords, n, alpha, beta, water=False) - mol_tension = smd.molecular_surface_tension(beta, gamma, phi, psi) + log = logger.new_logger(mol, mol.verbose) + t1 = log.init_timer() + if not smdobj._intermediates: + smdobj.build() + dm_cache = smdobj._intermediates.get('dm', None) + if dm_cache is not None and cupy.linalg.norm(dm_cache - dm) < 1e-10: + pass else: - grad_tension = atomic_surface_tension(symbols, coords, n, alpha, beta, water=True) - atm_tension = smd.atomic_surface_tension(symbols, coords, n, alpha, beta, water=True) - mol_tension = 0.0 - - # generate surface for calculating SASA - rad = radii.VDW + 0.4/radii.BOHR - surface = pcm.gen_surface(mol, ng=smdobj.sasa_ng, rad=rad) - _, grad_area = pcm_grad.get_dF_dA(surface) - area = surface['area'] - gridslice = surface['gslice_by_atom'] - SASA = cupy.asarray([cupy.sum(area[p0:p1], axis=0) for p0,p1, in gridslice]).get() - grad_SASA = cupy.asarray([cupy.sum(grad_area[p0:p1], axis=0) for p0,p1, in gridslice]).get() - SASA *= radii.BOHR**2 - grad_SASA *= radii.BOHR**2 - mol_cds = mol_tension * np.sum(grad_SASA, axis=0) / 1000 - grad_tension *= radii.BOHR - atm_cds = np.einsum('i,ijx->jx', SASA, grad_tension.get()) / 1000 - atm_cds+= np.einsum('i,ijx->jx', atm_tension.get(), grad_SASA) / 1000 - return (mol_cds + atm_cds)/hartree2kcal # hartree + smdobj._get_vind(dm) + + gridslice = smdobj.surface['gslice_by_atom'] + v_grids = smdobj._intermediates['v_grids'] + A = smdobj._intermediates['A'] + D = smdobj._intermediates['D'] + S = smdobj._intermediates['S'] + K = smdobj._intermediates['K'] + q = smdobj._intermediates['q'] + + vK_1 = cupy.linalg.solve(K.T, v_grids) + + dF, dA = pcm_grad.get_dF_dA(smdobj.surface) + + with_D = smdobj.method.upper() in ['IEF-PCM', 'IEFPCM', 'SS(V)PE', 'SMD'] + dD, dS, dSii = pcm_grad.get_dD_dS(smdobj.surface, dF, with_D=with_D, with_S=True) + DA = D*A + + epsilon = smdobj.eps + de = cupy.zeros([smdobj.mol.natm,3]) + + # same as IEF-PCM + dD = dD.transpose([2,0,1]) + dS = dS.transpose([2,0,1]) + dSii = dSii.transpose([2,0,1]) + dA = dA.transpose([2,0,1]) + def contract_bra(a, B, c): + ''' i,xij,j->jx ''' + tmp = contract('i,xij->xj', a, B) + return (tmp * c).T + + def contract_ket(a, B, c): + ''' i,xij,j->ix ''' + tmp = B.dot(c) + return (a*tmp).T + + # IEF-PCM and SS(V)PE formally are the same in gradient calculation + # dR = f_eps/(2*pi) * (dD*A + D*dA), + # dK = dS - f_eps/(2*pi) * (dD*A*S + D*dA*S + D*A*dS) + f_epsilon = (epsilon - 1.0)/(epsilon + 1.0) + fac = f_epsilon/(2.0*np.pi) + + Av = A*v_grids + de_dR = 0.5*fac * contract_ket(vK_1, dD, Av) + de_dR -= 0.5*fac * contract_bra(vK_1, dD, Av) + de_dR = cupy.asarray([cupy.sum(de_dR[p0:p1], axis=0) for p0,p1 in gridslice]) + + vK_1_D = vK_1.dot(D) + vK_1_Dv = vK_1_D * v_grids + de_dR += 0.5*fac * contract('j,xjn->nx', vK_1_Dv, dA) + + de_dS0 = 0.5*contract_ket(vK_1, dS, q) + de_dS0 -= 0.5*contract_bra(vK_1, dS, q) + de_dS0 = cupy.asarray([cupy.sum(de_dS0[p0:p1], axis=0) for p0,p1 in gridslice]) + + vK_1_q = vK_1 * q + de_dS0 += 0.5*contract('i,xin->nx', vK_1_q, dSii) + + vK_1_DA = cupy.dot(vK_1, DA) + de_dS1 = 0.5*contract_ket(vK_1_DA, dS, q) + de_dS1 -= 0.5*contract_bra(vK_1_DA, dS, q) + de_dS1 = cupy.asarray([cupy.sum(de_dS1[p0:p1], axis=0) for p0,p1 in gridslice]) + + vK_1_DAq = vK_1_DA*q + de_dS1 += 0.5*contract('j,xjn->nx', vK_1_DAq, dSii) + + Sq = cupy.dot(S,q) + ASq = A*Sq + de_dD = 0.5*contract_ket(vK_1, dD, ASq) + de_dD -= 0.5*contract_bra(vK_1, dD, ASq) + de_dD = cupy.asarray([cupy.sum(de_dD[p0:p1], axis=0) for p0,p1 in gridslice]) + + vK_1_D = cupy.dot(vK_1, D) + de_dA = 0.5*contract('j,xjn->nx', vK_1_D*Sq, dA) # 0.5*cupy.einsum('j,xjn,j->nx', vK_1_D, dA, Sq) + + de_dK = de_dS0 - fac * (de_dD + de_dA + de_dS1) + de += de_dR - de_dK + + t1 = log.timer_debug1('grad solver', *t1) + return de.get() def make_grad_object(grad_method): - '''For grad_method in vacuum, add nuclear gradients of solvent pcmobj''' + '''For grad_method in vacuum, add nuclear gradients of solvent smdobj''' if grad_method.base.with_solvent.frozen: raise RuntimeError('Frozen solvent model is not avialbe for energy gradients') @@ -233,6 +137,8 @@ def make_grad_object(grad_method): (WithSolventGrad, grad_method.__class__), name) class WithSolventGrad: + from gpu4pyscf.lib.utils import to_gpu, device + _keys = {'de_solvent', 'de_solute'} def __init__(self, grad_method): @@ -248,6 +154,11 @@ def undo_solvent(self): del obj.de_solute return obj + def to_cpu(self): + from pyscf.solvent.grad import smd # type: ignore + grad_method = self.undo_solvent().to_cpu() + return smd.make_grad_object(grad_method) + def kernel(self, *args, dm=None, atmlst=None, **kwargs): dm = kwargs.pop('dm', None) if dm is None: @@ -256,7 +167,7 @@ def kernel(self, *args, dm=None, atmlst=None, **kwargs): dm = dm[0] + dm[1] self.de_solute = super().kernel(*args, **kwargs) self.de_solvent = pcm_grad.grad_qv(self.base.with_solvent, dm) - self.de_solvent+= pcm_grad.grad_solver(self.base.with_solvent, dm) + self.de_solvent+= grad_solver(self.base.with_solvent, dm) self.de_solvent+= pcm_grad.grad_nuc(self.base.with_solvent, dm) self.de_cds = get_cds(self.base.with_solvent) self.de = self.de_solute + self.de_solvent + self.de_cds diff --git a/gpu4pyscf/solvent/grad/smd_experiment.py b/gpu4pyscf/solvent/grad/smd_experiment.py new file mode 100644 index 00000000..7836a005 --- /dev/null +++ b/gpu4pyscf/solvent/grad/smd_experiment.py @@ -0,0 +1,218 @@ +# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program 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 General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +''' +Gradient of SMD solvent model (for experiment and education) +''' + +import numpy as np +import cupy +from pyscf.data import radii +from gpu4pyscf.solvent import pcm +from gpu4pyscf.solvent import smd_experiment as smd +from gpu4pyscf.solvent.grad import pcm as pcm_grad + +def grad_switch_function(R, r, dr): + if R < r + dr: + return -np.exp(dr/(R-dr-r)) * dr / (R-dr-r)**2 + else: + return 0.0 + +from gpu4pyscf.solvent.smd import ( + sigma_water, sigma_n, sigma_alpha, sigma_beta, r_zz, switch_function, + hartree2kcal) + +def atomic_surface_tension(symbols, coords, n, alpha, beta, water=True): + ''' + - list of atomic symbols + - atomic coordinates in Anstrong + - solvent descriptors: n, alpha, beta + ''' + + def get_bond_tension(bond): + if water: + return sigma_water.get(bond, 0.0) + t = 0.0 + t += sigma_n.get(bond, 0.0) * n + t += sigma_alpha.get(bond, 0.0) * alpha + t += sigma_beta.get(bond, 0.0) * beta + return t + + def get_atom_tension(sym_i): + if water: + return sigma_water.get(sym_i, 0.0) + t = 0.0 + t += sigma_n.get(sym_i, 0.0) * n + t += sigma_alpha.get(sym_i, 0.0) * alpha + t += sigma_beta.get(sym_i, 0.0) * beta + return t + natm = coords.shape[0] + ri_rj = coords[:,None,:] - coords[None,:,:] + rij = np.sum(ri_rj**2, axis=2)**0.5 + np.fill_diagonal(rij, 1) + drij = ri_rj / np.expand_dims(rij, axis=-1) + tensions = [] + for i, sym_i in enumerate(symbols): + if sym_i not in ['H', 'C', 'N', 'O', 'F', 'Si', 'S', 'Cl', 'Br']: + tensions.append(np.zeros([natm,3])) + continue + + tension = np.zeros([natm,3]) + if sym_i in ['F', 'Si', 'S', 'Cl', 'Br']: + tensions.append(tension) + continue + + if sym_i == 'H': + dt_HC = np.zeros([natm,3]) + dt_HO = np.zeros([natm,3]) + for j, sym_j in enumerate(symbols): + if sym_j == 'C': + r, dr = r_zz.get(('H','C'), (0.0, 0.0)) + dt_drij = grad_switch_function(rij[i,j], r, dr) * drij[i,j] + dt_HC[i] += dt_drij + dt_HC[j] -= dt_drij + if sym_j == 'O': + r, dr = r_zz.get(('H','O'), (0.0, 0.0)) + dt_drij = grad_switch_function(rij[i,j], r, dr) * drij[i,j] + dt_HO[i] += dt_drij + dt_HO[j] -= dt_drij + sig_HC = get_bond_tension(('H','C')) + sig_HO = get_bond_tension(('H','O')) + tension += sig_HC * dt_HC + sig_HO * dt_HO + tensions.append(tension) + continue + + if sym_i == 'C': + dt_CC = np.zeros([natm,3]) + dt_CN = np.zeros([natm,3]) + t_CN = 0.0 + for j, sym_j in enumerate(symbols): + if sym_j == 'C' and i != j: + r, dr = r_zz.get(('C', 'C'), (0.0, 0.0)) + dt_drij = grad_switch_function(rij[i,j], r, dr) * drij[i,j] + dt_CC[i] += dt_drij + dt_CC[j] -= dt_drij + if sym_j == 'N': + r, dr = r_zz.get(('C', 'N'), (0.0, 0.0)) + t_CN += switch_function(rij[i,j], r, dr) + dt_drij = grad_switch_function(rij[i,j], r, dr) * drij[i,j] + dt_CN[i] += dt_drij + dt_CN[j] -= dt_drij + sig_CC = get_bond_tension(('C','C')) + sig_CN = get_bond_tension(('C','N')) + tension += sig_CC * dt_CC + sig_CN * (2 * t_CN * dt_CN) + tensions.append(tension) + continue + + if sym_i == 'N': + t_NC = 0.0 + dt_NC = np.zeros([natm,3]) + dt_NC3 = np.zeros([natm,3]) + for j, sym_j in enumerate(symbols): + if sym_j == 'C': + r, dr = r_zz.get(('N','C'), (0.0, 0.0)) + tk = 0.0 + dtk = np.zeros([natm,3]) + for k, sym_k in enumerate(symbols): + if k != i and k != j: + rjk, drjk = r_zz.get(('C', sym_k), (0.0, 0.0)) + tk += switch_function(rij[j,k], rjk, drjk) + dtk_rjk = grad_switch_function(rij[j,k], rjk, drjk) * drij[j,k] + dtk[j] += dtk_rjk + dtk[k] -= dtk_rjk + + dt_drij = grad_switch_function(rij[i,j], r, dr) * drij[i,j] * tk**2 + dt_NC[i] += dt_drij + dt_NC[j] -= dt_drij + + t = switch_function(rij[i,j], r, dr) + dt_NC += t * (2 * tk * dtk) + t_NC += t * tk**2 + + r, dr = r_zz.get(('N','C3'), (0.0, 0.0)) + dt_drij = grad_switch_function(rij[i,j], r, dr) * drij[i,j] + dt_NC3[i] += dt_drij + dt_NC3[j] -= dt_drij + sig_NC = get_bond_tension(('N','C')) + sig_NC3= get_bond_tension(('N','C3')) + tension += sig_NC * (1.3 * t_NC**0.3 * dt_NC) + sig_NC3 * dt_NC3 + tensions.append(tension) + continue + + if sym_i == 'O': + dt_OC = np.zeros([natm,3]) + dt_ON = np.zeros([natm,3]) + dt_OO = np.zeros([natm,3]) + dt_OP = np.zeros([natm,3]) + for j, sym_j in enumerate(symbols): + if sym_j == 'C': + r, dr = r_zz.get(('O','C'), (0.0, 0.0)) + dt_drij = grad_switch_function(rij[i,j], r, dr) * drij[i,j] + dt_OC[i] += dt_drij + dt_OC[j] -= dt_drij + if sym_j == 'N': + r, dr = r_zz.get(('O','N'), (0.0, 0.0)) + dt_drij = grad_switch_function(rij[i,j], r, dr) * drij[i,j] + dt_ON[i] += dt_drij + dt_ON[j] -= dt_drij + if sym_j == 'O' and j != i: + r, dr = r_zz.get(('O','O'), (0.0, 0.0)) + dt_drij = grad_switch_function(rij[i,j], r, dr) * drij[i,j] + dt_OO[i] += dt_drij + dt_OO[j] -= dt_drij + if sym_j == 'P': + r, dr = r_zz.get(('O','P'), (0.0, 0.0)) + dt_drij = grad_switch_function(rij[i,j], r, dr) * drij[i,j] + dt_OP[i] += dt_drij + dt_OP[j] -= dt_drij + sig_OC = get_bond_tension(('O','C')) + sig_ON = get_bond_tension(('O','N')) + sig_OO = get_bond_tension(('O','O')) + sig_OP = get_bond_tension(('O','P')) + tension += sig_OC * dt_OC + sig_ON * dt_ON + sig_OO * dt_OO + sig_OP * dt_OP + tensions.append(tension) + continue + + return cupy.asarray(tensions) + +def get_cds(smdobj): + mol = smdobj.mol + n, _, alpha, beta, gamma, _, phi, psi = smdobj.solvent_descriptors + symbols = [mol.atom_symbol(ia) for ia in range(mol.natm)] + coords = mol.atom_coords(unit='A') + if smdobj._solvent.lower() != 'water': + grad_tension = atomic_surface_tension(symbols, coords, n, alpha, beta, water=False) + atm_tension = smd.atomic_surface_tension(symbols, coords, n, alpha, beta, water=False) + mol_tension = smd.molecular_surface_tension(beta, gamma, phi, psi) + else: + grad_tension = atomic_surface_tension(symbols, coords, n, alpha, beta, water=True) + atm_tension = smd.atomic_surface_tension(symbols, coords, n, alpha, beta, water=True) + mol_tension = 0.0 + + # generate surface for calculating SASA + rad = radii.VDW + 0.4/radii.BOHR + surface = pcm.gen_surface(mol, ng=smdobj.sasa_ng, rad=rad) + _, grad_area = pcm_grad.get_dF_dA(surface) + area = surface['area'] + gridslice = surface['gslice_by_atom'] + SASA = cupy.asarray([cupy.sum(area[p0:p1], axis=0) for p0,p1, in gridslice]).get() + grad_SASA = cupy.asarray([cupy.sum(grad_area[p0:p1], axis=0) for p0,p1, in gridslice]).get() + SASA *= radii.BOHR**2 + grad_SASA *= radii.BOHR**2 + mol_cds = mol_tension * np.sum(grad_SASA, axis=0) / 1000 + grad_tension *= radii.BOHR + atm_cds = np.einsum('i,ijx->jx', SASA, grad_tension.get()) / 1000 + atm_cds+= np.einsum('i,ijx->jx', atm_tension.get(), grad_SASA) / 1000 + return (mol_cds + atm_cds)/hartree2kcal # hartree \ No newline at end of file diff --git a/gpu4pyscf/solvent/hessian/pcm.py b/gpu4pyscf/solvent/hessian/pcm.py index 840a8090..c6473fcf 100644 --- a/gpu4pyscf/solvent/hessian/pcm.py +++ b/gpu4pyscf/solvent/hessian/pcm.py @@ -14,7 +14,7 @@ # along with this program. If not, see . ''' -Gradient of PCM family solvent model +Hessian of PCM family solvent model ''' # pylint: disable=C0103 @@ -22,14 +22,12 @@ import cupy from pyscf import lib, gto from gpu4pyscf import scf -from gpu4pyscf.solvent.pcm import PI, switch_h -from gpu4pyscf.solvent.grad.pcm import grad_switch_h, get_dF_dA, get_dD_dS, grad_qv, grad_solver, grad_nuc +from gpu4pyscf.solvent.pcm import PI +from gpu4pyscf.solvent.grad.pcm import grad_qv, grad_solver, grad_nuc from gpu4pyscf.df import int3c2e from gpu4pyscf.lib.cupy_helper import contract from gpu4pyscf.lib import logger -libdft = lib.load_library('libdft') - def hess_nuc(pcmobj): if not pcmobj._intermediates: pcmobj.build() @@ -142,11 +140,9 @@ def pcm_grad_scanner(mol): dv = numpy.zeros_like(coords) dv[ia,ix] = eps mol.set_geom_(coords + dv, unit='Bohr') - mol.build() g0 = pcm_grad_scanner(mol) mol.set_geom_(coords - dv, unit='Bohr') - mol.build() g1 = pcm_grad_scanner(mol) de[ia,:,ix] = (g0 - g1)/2.0/eps t1 = log.timer_debug1('solvent energy', *t1) @@ -181,11 +177,9 @@ def pcm_vmat_scanner(mol): dv = numpy.zeros_like(coords) dv[ia,ix] = eps mol.set_geom_(coords + dv, unit='Bohr') - mol.build() vmat0 = pcm_vmat_scanner(mol) mol.set_geom_(coords - dv, unit='Bohr') - mol.build() vmat1 = pcm_vmat_scanner(mol) grad_vmat = (vmat0 - vmat1)/2.0/eps @@ -236,6 +230,8 @@ def make_hess_object(hess_method): (WithSolventHess, hess_method.__class__), name) class WithSolventHess: + from gpu4pyscf.lib.utils import to_gpu, device + _keys = {'de_solvent', 'de_solute'} def __init__(self, hess_method): @@ -243,6 +239,19 @@ def __init__(self, hess_method): self.de_solvent = None self.de_solute = None + def undo_solvent(self): + cls = self.__class__ + name_mixin = self.base.with_solvent.__class__.__name__ + obj = lib.view(self, lib.drop_class(cls, WithSolventHess, name_mixin)) + del obj.de_solvent + del obj.de_solute + return obj + + def to_cpu(self): + from pyscf.solvent.hessian import pcm # type: ignore + hess_method = self.undo_solvent().to_cpu() + return pcm.make_hess_object(hess_method) + def kernel(self, *args, dm=None, atmlst=None, **kwargs): dm = kwargs.pop('dm', None) if dm is None: diff --git a/gpu4pyscf/solvent/hessian/smd.py b/gpu4pyscf/solvent/hessian/smd.py index a97e0d14..3e8c5238 100644 --- a/gpu4pyscf/solvent/hessian/smd.py +++ b/gpu4pyscf/solvent/hessian/smd.py @@ -14,216 +14,26 @@ # along with this program. If not, see . ''' -Gradient of PCM family solvent model +Hessian SMD solvent model ''' # pylint: disable=C0103 import numpy as np -import cupy from pyscf import lib - +from gpu4pyscf import scf +from gpu4pyscf.lib import logger from gpu4pyscf.solvent import smd from gpu4pyscf.solvent.grad import smd as smd_grad +from gpu4pyscf.solvent.grad import pcm as pcm_grad from gpu4pyscf.solvent.hessian import pcm as pcm_hess -from gpu4pyscf.solvent.smd import ( - sigma_water, sigma_n, sigma_alpha, sigma_beta, r_zz, swtich_function, - hartree2kcal) -from gpu4pyscf import scf -from gpu4pyscf.lib.cupy_helper import contract -from gpu4pyscf.lib import logger - -def hess_swtich_function(R, r, dr): - if R < r + dr: - dist = (R[0]**2 + R[1]**2 + R[2]**2)**.5 - hess = [ - [R[1]**2+R[2]**2, -R[0]*R[1], -R[0]*R[2]], - [-R[0]*R[1], R[0]**2+R[2]**2, -R[1]*R[2]], - [-R[0]*R[2], -R[1]*R[2], R[0]**2+R[1]**2]] - return np.asarray(hess)/dist - else: - return np.zeros([3,3]) - -def atomic_surface_tension(symbols, coords, n, alpha, beta, water=True): - ''' - TODO: debug later - - list of atomic symbols - - atomic coordinates in Anstrong - - solvent descriptors: n, alpha, beta - ''' - - def get_bond_tension(bond): - if water: - return sigma_water.get(bond, 0.0) - t = 0.0 - t += sigma_n.get(bond, 0.0) * n - t += sigma_alpha.get(bond, 0.0) * alpha - t += sigma_beta.get(bond, 0.0) * beta - return t - - def get_atom_tension(sym_i): - if water: - return sigma_water.get(sym_i, 0.0) - t = 0.0 - t += sigma_n.get(sym_i, 0.0) * n - t += sigma_alpha.get(sym_i, 0.0) * alpha - t += sigma_beta.get(sym_i, 0.0) * beta - return t - natm = coords.shape[0] - ri_rj = coords[:,None,:] - coords[None,:,:] - rij = np.sum(ri_rj**2, axis=2)**0.5 - np.fill_diagonal(rij, 1) - #drij = ri_rj / np.expand_dims(rij, axis=-1) - tensions = [] - for i, sym_i in enumerate(symbols): - if sym_i not in ['H', 'C', 'N', 'O', 'F', 'Si', 'S', 'Cl', 'Br']: - tensions.append(np.zeros([natm,3])) - continue - - tension = np.zeros([natm,3]) - if sym_i in ['F', 'Si', 'S', 'Cl', 'Br']: - tensions.append(tension) - continue - - if sym_i == 'H': - dt_HC = np.zeros([natm,natm,3,3]) - dt_HO = np.zeros([natm,natm,3,3]) - for j, sym_j in enumerate(symbols): - if sym_j == 'C': - r, dr = r_zz.get(('H','C'), (0.0, 0.0)) - dt_drij = hess_swtich_function(coords[i]-coords[j], r, dr) - dt_HC[i,i] += dt_drij - dt_HC[i,j] -= dt_drij - dt_HC[j,i] -= dt_drij - dt_HC[j,j] += dt_drij - if sym_j == 'O': - r, dr = r_zz.get(('H','O'), (0.0, 0.0)) - dt_drij = hess_swtich_function(coords[i]-coords[j], r, dr) - dt_HO[i,i] += dt_drij - dt_HO[i,j] -= dt_drij - dt_HO[j,i] -= dt_drij - dt_HO[j,j] += dt_drij - sig_HC = get_bond_tension(('H','C')) - sig_HO = get_bond_tension(('H','O')) - tension += sig_HC * dt_HC + sig_HO * dt_HO - tensions.append(tension) - - if sym_i == 'C': - dt_CN = np.zeros([natm,3]) - d2t_CC = np.zeros([natm,natm,3,3]) - d2t_CN = np.zeros([natm,natm,3,3]) - t_CN = 0.0 - for j, sym_j in enumerate(symbols): - if sym_j == 'C' and i != j: - r, dr = r_zz.get(('C', 'C'), (0.0, 0.0)) - d2t_drij = hess_swtich_function(coords[i]-coords[j], r, dr) - d2t_CC[i,i] += d2t_drij - d2t_CC[i,j] -= d2t_drij - d2t_CC[j,i] -= d2t_drij - d2t_CC[j,j] += d2t_drij - if sym_j == 'N': - r, dr = r_zz.get(('C', 'N'), (0.0, 0.0)) - t_CN += swtich_function(rij[i,j], r, dr) - dt_drij = smd_grad.grad_switch_function(coords[i]-coords[j], r, dr) - dt_CN[i] += dt_drij - dt_CN[j] -= dt_drij - d2t_drij = hess_swtich_function(coords[i]-coords[j], r, dr) - d2t_CN[i,i] += d2t_drij - d2t_CN[i,j] -= d2t_drij - d2t_CN[j,i] -= d2t_drij - d2t_CN[j,j] += d2t_drij - sig_CC = get_bond_tension(('C','C')) - sig_CN = get_bond_tension(('C','N')) - tension += sig_CC * d2t_CC + sig_CN * (2 * t_CN * d2t_CN + 2 * np.einsum('i,j->ij', dt_drij[i], dt_drij[j])) - tensions.append(tension) - - if sym_i == 'N': - t_NC = 0.0 - dt_NC = np.zeros([natm,natm,3,3]) - dt_NC3 = np.zeros([natm,natm,3,3]) - for j, sym_j in enumerate(symbols): - if sym_j == 'C': - r, dr = r_zz.get(('N','C'), (0.0, 0.0)) - tk = 0.0 - dtk = np.zeros([natm,natm,3,3]) - for k, sym_k in enumerate(symbols): - if k != i and k != j: - rjk, drjk = r_zz.get(('C', sym_k), (0.0, 0.0)) - tk += swtich_function(rij[j,k], rjk, drjk) - dtk_rjk = hess_swtich_function(coords[j]-coords[k], rjk, drjk) - dtk[j,j] += dtk_rjk - dtk[j,k] -= dtk_rjk - dtk[k,j] -= dtk_rjk - dtk[k,k] += dtk_rjk - dt_drij = hess_swtich_function(coords[i]-coords[j], r, dr) * tk**2 - dt_NC[i,i] += dt_drij - dt_NC[i,j] -= dt_drij - dt_NC[j,i] -= dt_drij - dt_NC[j,j] += dt_drij - t = swtich_function(coords[i]-coords[j], r, dr) - dt_NC += t * (2 * tk * dtk) - t_NC += t * tk**2 - - r, dr = r_zz.get(('N','C3'), (0.0, 0.0)) - dt_drij = hess_swtich_function(coords[i]-coords[j], r, dr) - dt_NC3[i,i] += dt_drij - dt_NC3[i,j] -= dt_drij - dt_NC3[j,i] -= dt_drij - dt_NC3[j,j] += dt_drij - sig_NC = get_bond_tension(('N','C')) - sig_NC3= get_bond_tension(('N','C3')) - tension += sig_NC * (1.3 * t_NC**0.3 * dt_NC) + sig_NC3 * dt_NC3 - tensions.append(tension) - - if sym_i == 'O': - dt_OC = np.zeros([natm,natm,3,3]) - dt_ON = np.zeros([natm,natm,3,3]) - dt_OO = np.zeros([natm,natm,3,3]) - dt_OP = np.zeros([natm,natm,3,3]) - for j, sym_j in enumerate(symbols): - if sym_j == 'C': - r, dr = r_zz.get(('O','C'), (0.0, 0.0)) - dt_drij = hess_swtich_function(coords[i]-coords[j], r, dr) - dt_OC[i,i] += dt_drij - dt_OC[i,j] -= dt_drij - dt_OC[j,i] -= dt_drij - dt_OC[j,j] += dt_drij - if sym_j == 'N': - r, dr = r_zz.get(('O','N'), (0.0, 0.0)) - dt_drij = hess_swtich_function(coords[i]-coords[j], r, dr) - dt_ON[i,i] += dt_drij - dt_ON[i,j] -= dt_drij - dt_ON[j,i] -= dt_drij - dt_ON[j,j] += dt_drij - if sym_j == 'O' and j != i: - r, dr = r_zz.get(('O','O'), (0.0, 0.0)) - dt_drij = hess_swtich_function(coords[i]-coords[j], r, dr) - dt_OO[i,i] += dt_drij - dt_OO[i,j] -= dt_drij - dt_OO[j,i] -= dt_drij - dt_OO[j,j] += dt_drij - if sym_j == 'P': - r, dr = r_zz.get(('O','P'), (0.0, 0.0)) - dt_drij = hess_swtich_function(coords[i]-coords[j], r, dr) - dt_OP[i,i] += dt_drij - dt_OP[i,j] -= dt_drij - dt_OP[j,i] -= dt_drij - dt_OP[j,j] += dt_drij - sig_OC = get_bond_tension(('O','C')) - sig_ON = get_bond_tension(('O','N')) - sig_OO = get_bond_tension(('O','O')) - sig_OP = get_bond_tension(('O','P')) - tension += sig_OC * dt_OC + sig_ON * dt_ON + sig_OO * dt_OO + sig_OP * dt_OP - tensions.append(tension) - return cupy.asarray(tensions) - def get_cds(smdobj): mol = smdobj.mol solvent = smdobj.solvent def smd_grad_scanner(mol): smdobj_tmp = smd.SMD(mol) smdobj_tmp.solvent = solvent - return smd_grad.get_cds(smdobj_tmp) + return smd.get_cds_legacy(smdobj_tmp)[1] log = logger.new_logger(mol, mol.verbose) t1 = log.init_timer() @@ -250,6 +60,45 @@ def smd_grad_scanner(mol): t1 = log.timer_debug1('solvent energy', *t1) return hess_cds # hartree + +def hess_elec(smdobj, dm, verbose=None): + ''' + slow version with finite difference + TODO: use analytical hess_nuc + ''' + log = logger.new_logger(smdobj, verbose) + t1 = log.init_timer() + pmol = smdobj.mol.copy() + mol = pmol.copy() + coords = mol.atom_coords(unit='Bohr') + + def pcm_grad_scanner(mol): + # TODO: use more analytical forms + smdobj.reset(mol) + e, v = smdobj._get_vind(dm) + #return grad_elec(smdobj, dm) + grad = pcm_grad.grad_nuc(smdobj, dm) + grad+= smd_grad.grad_solver(smdobj, dm) + grad+= pcm_grad.grad_qv(smdobj, dm) + return grad + + mol.verbose = 0 + de = np.zeros([mol.natm, mol.natm, 3, 3]) + eps = 1e-3 + for ia in range(mol.natm): + for ix in range(3): + dv = np.zeros_like(coords) + dv[ia,ix] = eps + mol.set_geom_(coords + dv, unit='Bohr') + g0 = pcm_grad_scanner(mol) + + mol.set_geom_(coords - dv, unit='Bohr') + g1 = pcm_grad_scanner(mol) + de[ia,:,ix] = (g0 - g1)/2.0/eps + t1 = log.timer_debug1('solvent energy', *t1) + smdobj.reset(pmol) + return de + def make_hess_object(hess_method): '''For hess_method in vacuum, add nuclear Hessian of solvent smdobj''' if hess_method.base.with_solvent.frozen: @@ -261,6 +110,8 @@ def make_hess_object(hess_method): (WithSolventHess, hess_method.__class__), name) class WithSolventHess: + from gpu4pyscf.lib.utils import to_gpu, device + _keys = {'de_solvent', 'de_solute'} def __init__(self, hess_method): @@ -268,6 +119,19 @@ def __init__(self, hess_method): self.de_solvent = None self.de_solute = None + def undo_solvent(self): + cls = self.__class__ + name_mixin = self.base.with_solvent.__class__.__name__ + obj = lib.view(self, lib.drop_class(cls, WithSolventHess, name_mixin)) + del obj.de_solvent + del obj.de_solute + return obj + + def to_cpu(self): + from pyscf.solvent.hessian import smd # type: ignore + hess_method = self.undo_solvent().to_cpu() + return smd.make_hess_object(hess_method) + def kernel(self, *args, dm=None, atmlst=None, **kwargs): dm = kwargs.pop('dm', None) if dm is None: diff --git a/gpu4pyscf/solvent/hessian/smd_experiment.py b/gpu4pyscf/solvent/hessian/smd_experiment.py new file mode 100644 index 00000000..f90cfe38 --- /dev/null +++ b/gpu4pyscf/solvent/hessian/smd_experiment.py @@ -0,0 +1,193 @@ +# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program 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 General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +''' +Hessian SMD solvent model (for experiment and education) +''' + +import cupy +import numpy as np +from gpu4pyscf.solvent.grad import smd as smd_grad +from gpu4pyscf.solvent.smd_experiment import ( + sigma_water, sigma_n, sigma_alpha, sigma_beta, r_zz, switch_function) + +def hess_switch_function(R, r, dr): + if R < r + dr: + dist = (R[0]**2 + R[1]**2 + R[2]**2)**.5 + hess = [ + [R[1]**2+R[2]**2, -R[0]*R[1], -R[0]*R[2]], + [-R[0]*R[1], R[0]**2+R[2]**2, -R[1]*R[2]], + [-R[0]*R[2], -R[1]*R[2], R[0]**2+R[1]**2]] + return np.asarray(hess)/dist + else: + return np.zeros([3,3]) + +def atomic_surface_tension(symbols, coords, n, alpha, beta, water=True): + def get_bond_tension(bond): + if water: + return sigma_water.get(bond, 0.0) + t = 0.0 + t += sigma_n.get(bond, 0.0) * n + t += sigma_alpha.get(bond, 0.0) * alpha + t += sigma_beta.get(bond, 0.0) * beta + return t + + natm = coords.shape[0] + ri_rj = coords[:,None,:] - coords[None,:,:] + rij = np.sum(ri_rj**2, axis=2)**0.5 + np.fill_diagonal(rij, 1) + #drij = ri_rj / np.expand_dims(rij, axis=-1) + tensions = [] + for i, sym_i in enumerate(symbols): + if sym_i not in ['H', 'C', 'N', 'O', 'F', 'Si', 'S', 'Cl', 'Br']: + tensions.append(np.zeros([natm,3])) + continue + + tension = np.zeros([natm,3]) + if sym_i in ['F', 'Si', 'S', 'Cl', 'Br']: + tensions.append(tension) + continue + + if sym_i == 'H': + dt_HC = np.zeros([natm,natm,3,3]) + dt_HO = np.zeros([natm,natm,3,3]) + for j, sym_j in enumerate(symbols): + if sym_j == 'C': + r, dr = r_zz.get(('H','C'), (0.0, 0.0)) + dt_drij = hess_switch_function(coords[i]-coords[j], r, dr) + dt_HC[i,i] += dt_drij + dt_HC[i,j] -= dt_drij + dt_HC[j,i] -= dt_drij + dt_HC[j,j] += dt_drij + if sym_j == 'O': + r, dr = r_zz.get(('H','O'), (0.0, 0.0)) + dt_drij = hess_switch_function(coords[i]-coords[j], r, dr) + dt_HO[i,i] += dt_drij + dt_HO[i,j] -= dt_drij + dt_HO[j,i] -= dt_drij + dt_HO[j,j] += dt_drij + sig_HC = get_bond_tension(('H','C')) + sig_HO = get_bond_tension(('H','O')) + tension += sig_HC * dt_HC + sig_HO * dt_HO + tensions.append(tension) + + if sym_i == 'C': + dt_CN = np.zeros([natm,3]) + d2t_CC = np.zeros([natm,natm,3,3]) + d2t_CN = np.zeros([natm,natm,3,3]) + t_CN = 0.0 + for j, sym_j in enumerate(symbols): + if sym_j == 'C' and i != j: + r, dr = r_zz.get(('C', 'C'), (0.0, 0.0)) + d2t_drij = hess_switch_function(coords[i]-coords[j], r, dr) + d2t_CC[i,i] += d2t_drij + d2t_CC[i,j] -= d2t_drij + d2t_CC[j,i] -= d2t_drij + d2t_CC[j,j] += d2t_drij + if sym_j == 'N': + r, dr = r_zz.get(('C', 'N'), (0.0, 0.0)) + t_CN += switch_function(rij[i,j], r, dr) + dt_drij = smd_grad.grad_switch_function(coords[i]-coords[j], r, dr) + dt_CN[i] += dt_drij + dt_CN[j] -= dt_drij + d2t_drij = hess_switch_function(coords[i]-coords[j], r, dr) + d2t_CN[i,i] += d2t_drij + d2t_CN[i,j] -= d2t_drij + d2t_CN[j,i] -= d2t_drij + d2t_CN[j,j] += d2t_drij + sig_CC = get_bond_tension(('C','C')) + sig_CN = get_bond_tension(('C','N')) + tension += sig_CC * d2t_CC + sig_CN * (2 * t_CN * d2t_CN + 2 * np.einsum('i,j->ij', dt_drij[i], dt_drij[j])) + tensions.append(tension) + + if sym_i == 'N': + t_NC = 0.0 + dt_NC = np.zeros([natm,natm,3,3]) + dt_NC3 = np.zeros([natm,natm,3,3]) + for j, sym_j in enumerate(symbols): + if sym_j == 'C': + r, dr = r_zz.get(('N','C'), (0.0, 0.0)) + tk = 0.0 + dtk = np.zeros([natm,natm,3,3]) + for k, sym_k in enumerate(symbols): + if k != i and k != j: + rjk, drjk = r_zz.get(('C', sym_k), (0.0, 0.0)) + tk += switch_function(rij[j,k], rjk, drjk) + dtk_rjk = hess_switch_function(coords[j]-coords[k], rjk, drjk) + dtk[j,j] += dtk_rjk + dtk[j,k] -= dtk_rjk + dtk[k,j] -= dtk_rjk + dtk[k,k] += dtk_rjk + dt_drij = hess_switch_function(coords[i]-coords[j], r, dr) * tk**2 + dt_NC[i,i] += dt_drij + dt_NC[i,j] -= dt_drij + dt_NC[j,i] -= dt_drij + dt_NC[j,j] += dt_drij + t = switch_function(coords[i]-coords[j], r, dr) + dt_NC += t * (2 * tk * dtk) + t_NC += t * tk**2 + + r, dr = r_zz.get(('N','C3'), (0.0, 0.0)) + dt_drij = hess_switch_function(coords[i]-coords[j], r, dr) + dt_NC3[i,i] += dt_drij + dt_NC3[i,j] -= dt_drij + dt_NC3[j,i] -= dt_drij + dt_NC3[j,j] += dt_drij + sig_NC = get_bond_tension(('N','C')) + sig_NC3= get_bond_tension(('N','C3')) + tension += sig_NC * (1.3 * t_NC**0.3 * dt_NC) + sig_NC3 * dt_NC3 + tensions.append(tension) + + if sym_i == 'O': + dt_OC = np.zeros([natm,natm,3,3]) + dt_ON = np.zeros([natm,natm,3,3]) + dt_OO = np.zeros([natm,natm,3,3]) + dt_OP = np.zeros([natm,natm,3,3]) + for j, sym_j in enumerate(symbols): + if sym_j == 'C': + r, dr = r_zz.get(('O','C'), (0.0, 0.0)) + dt_drij = hess_switch_function(coords[i]-coords[j], r, dr) + dt_OC[i,i] += dt_drij + dt_OC[i,j] -= dt_drij + dt_OC[j,i] -= dt_drij + dt_OC[j,j] += dt_drij + if sym_j == 'N': + r, dr = r_zz.get(('O','N'), (0.0, 0.0)) + dt_drij = hess_switch_function(coords[i]-coords[j], r, dr) + dt_ON[i,i] += dt_drij + dt_ON[i,j] -= dt_drij + dt_ON[j,i] -= dt_drij + dt_ON[j,j] += dt_drij + if sym_j == 'O' and j != i: + r, dr = r_zz.get(('O','O'), (0.0, 0.0)) + dt_drij = hess_switch_function(coords[i]-coords[j], r, dr) + dt_OO[i,i] += dt_drij + dt_OO[i,j] -= dt_drij + dt_OO[j,i] -= dt_drij + dt_OO[j,j] += dt_drij + if sym_j == 'P': + r, dr = r_zz.get(('O','P'), (0.0, 0.0)) + dt_drij = hess_switch_function(coords[i]-coords[j], r, dr) + dt_OP[i,i] += dt_drij + dt_OP[i,j] -= dt_drij + dt_OP[j,i] -= dt_drij + dt_OP[j,j] += dt_drij + sig_OC = get_bond_tension(('O','C')) + sig_ON = get_bond_tension(('O','N')) + sig_OO = get_bond_tension(('O','O')) + sig_OP = get_bond_tension(('O','P')) + tension += sig_OC * dt_OC + sig_ON * dt_ON + sig_OO * dt_OO + sig_OP * dt_OP + tensions.append(tension) + return cupy.asarray(tensions) diff --git a/gpu4pyscf/solvent/pcm.py b/gpu4pyscf/solvent/pcm.py index 39bd17b5..0e7fdd97 100644 --- a/gpu4pyscf/solvent/pcm.py +++ b/gpu4pyscf/solvent/pcm.py @@ -22,10 +22,10 @@ import cupy import cupyx.scipy as scipy from pyscf import lib -from pyscf import gto, df +from pyscf import gto from pyscf.dft import gen_grid from pyscf.data import radii -from pyscf.solvent import ddcosmo#, _attach_solvent +from pyscf.solvent import ddcosmo from gpu4pyscf.solvent import _attach_solvent from gpu4pyscf.df import int3c2e from gpu4pyscf.lib import logger @@ -244,33 +244,53 @@ class PCM(lib.StreamObject): 'eps', 'grids', 'max_cycle', 'conv_tol', 'state_id', 'frozen', 'equilibrium_solvation', 'e', 'v', } - + from gpu4pyscf.lib.utils import to_gpu, device kernel = ddcosmo.DDCOSMO.kernel def __init__(self, mol): - ddcosmo.DDCOSMO.__init__(self, mol) + self.mol = mol + self.stdout = mol.stdout + self.verbose = mol.verbose + self.max_memory = mol.max_memory self.method = 'C-PCM' + self.vdw_scale = 1.2 # default value in qchem + self.surface = {} self.r_probe = 0.0 self.radii_table = None - self.surface = {} + self.atom_radii = None + self.lebedev_order = 29 self._intermediates = {} + self.eps = 78.3553 + + self.max_cycle = 20 + self.conv_tol = 1e-7 + self.state_id = 0 + + self.frozen = False + self.equilibrium_solvation = False + + self.e = None + self.v = None + self._dm = None def dump_flags(self, verbose=None): logger.info(self, '******** %s ********', self.__class__) logger.info(self, 'lebedev_order = %s (%d grids per sphere)', self.lebedev_order, gen_grid.LEBEDEV_ORDER[self.lebedev_order]) - logger.info(self, 'lmax = %s' , self.lmax) - logger.info(self, 'eta = %s' , self.eta) logger.info(self, 'eps = %s' , self.eps) logger.info(self, 'frozen = %s' , self.frozen) logger.info(self, 'equilibrium_solvation = %s', self.equilibrium_solvation) logger.debug2(self, 'radii_table %s', self.radii_table) if self.atom_radii: logger.info(self, 'User specified atomic radii %s', str(self.atom_radii)) - self.grids.dump_flags(verbose) return self + def to_cpu(self): + from gpu4pyscf.lib.utils import to_cpu + obj = to_cpu(self) + return obj.reset() + def build(self, ng=None): if self.radii_table is None: vdw_scale = self.vdw_scale @@ -293,7 +313,7 @@ def build(self, ng=None): f_epsilon = (epsilon - 1.0)/(epsilon + 1.0/2.0) K = S R = -f_epsilon * cupy.eye(K.shape[0]) - elif self.method.upper() in ['IEF-PCM', 'IEFPCM', 'SMD']: + elif self.method.upper() in ['IEF-PCM', 'IEFPCM']: f_epsilon = (epsilon - 1.0)/(epsilon + 1.0) DA = D*A DAS = cupy.dot(DA, S) @@ -406,7 +426,6 @@ def Hessian(self, hess_method): def reset(self, mol=None): if mol is not None: self.mol = mol - self.grids.reset(mol) self._intermediates = None self.surface = None self.intopt = None diff --git a/gpu4pyscf/solvent/smd.py b/gpu4pyscf/solvent/smd.py index 06672e98..61704635 100644 --- a/gpu4pyscf/solvent/smd.py +++ b/gpu4pyscf/solvent/smd.py @@ -18,14 +18,13 @@ ''' import numpy as np -import scipy import cupy -from pyscf import lib +from pyscf import lib, gto from pyscf.data import radii from pyscf.dft import gen_grid from gpu4pyscf.solvent import pcm, _attach_solvent from gpu4pyscf.lib import logger -from gpu4pyscf.lib.cupy_helper import dist_matrix +from gpu4pyscf.df import int3c2e @lib.with_doc(_attach_solvent._for_scf.__doc__) def smd_for_scf(mf, solvent_obj=None, dm=None): @@ -37,82 +36,7 @@ def smd_for_scf(mf, solvent_obj=None, dm=None): from gpu4pyscf import scf scf.hf.RHF.SMD = smd_for_scf scf.uhf.UHF.SMD = smd_for_scf -hartree2kcal = 627.5 -# see https://pubs.acs.org/doi/epdf/10.1021/jp810292n - -sigma_water = { - ('H'): 48.69, - ('C'): 129.74, - ('H','C'): -60.77, - ('C','C'): -72.95, - ('O','C'): 68.69, - ('N','C'): -48.22, - ('N','C3'): 84.10, - ('O','N'): 121.98, - ('F'): 38.18, - ('Cl'): 9.82, - ('Br'): -8.72, - ('S'): -9.10, - ('O','P'): 68.85} - -sigma_n = { - ('C'): 58.10, - ('H','C'): -36.37, - ('C','C'): -62.05, - ('O'): -17.56, - ('H','O'): -19.39, - ('O','C'): -15.70, - ('N'): 32.62, - ('C','N'): -99.76, - ('Cl'): -24.31, - ('Br'): -35.42, - ('S'): -33.17, - ('Si'): -18.04 -} - -sigma_alpha = { - ('C'): 48.10, - ('O'): 193.06, - ('O','C'): 95.99, - ('C','N'): 152.20, - ('N','C'): -41.00 -} - -sigma_beta = { - ('C'): 32.87, - ('O'): -43.79, - ('O','O'):-128.16, - ('O','N'):79.13 -} - -# Molecular surface tension (cal/mol*AA^-2) -sigma_gamma = 0.35 -sigma_phi2 = -4.19 -sigma_psi2 = -6.68 -sigma_beta2 = 0.0 -gamma0 = 1.0 - -# rzz and delta r_zz in AA -r_zz = { - ('H','C'): [1.55, 0.3], - ('H','O'): [1.55, 0.3], - ('C','H'): [1.55, 0.3], - ('C','C'): [1.84, 0.3], - ('C','N'): [1.84, 0.3], - ('C','O'): [1.84, 0.3], - ('C','F'): [1.84, 0.3], - ('C','P'): [2.2, 0.3], - ('C','S'): [2.2, 0.3], - ('C','Cl'): [2.1, 0.3], - ('C','Br'): [2.3, 0.3], - ('C','I'): [2.6, 0.3], - ('N','C'): [1.84, 0.3], - ('N','C3'): [1.225, 0.065], - ('O','C'): [1.33, 0.1], - ('O','N'): [1.5, 0.3], - ('O','O'): [1.8, 0.3], - ('O','P'): [2.1, 0.3] -} +hartree2kcal = 627.509451 # database from https://comp.chem.umn.edu/solvation/mnsddb.pdf # solvent name: [n, n25, alpha, beta, gamma, epsilon, phi, psi) @@ -327,175 +251,48 @@ def smd_radii(alpha): radii_table[53] = 2.74 return radii_table/radii.BOHR -def swtich_function(R, r, dr): - return np.exp(dr/(R-dr-r)) if R. + +''' +SMD solvent model (for experiment and education) +''' + +import numpy as np +import scipy +import cupy +from pyscf.data import radii +from gpu4pyscf.solvent.smd import hartree2kcal +from gpu4pyscf.solvent import pcm +from gpu4pyscf.lib import logger + +# see https://pubs.acs.org/doi/epdf/10.1021/jp810292n + +sigma_water = { + ('H'): 48.69, + ('C'): 129.74, + ('H','C'): -60.77, + ('C','C'): -72.95, + ('O','C'): 68.69, + ('N','C'): -48.22, + ('N','C3'): 84.10, + ('O','N'): 121.98, + ('F'): 38.18, + ('Cl'): 9.82, + ('Br'): -8.72, + ('S'): -9.10, + ('O','P'): 68.85} + +sigma_n = { + ('C'): 58.10, + ('H','C'): -36.37, + ('C','C'): -62.05, + ('O'): -17.56, + ('H','O'): -19.39, + ('O','C'): -15.70, + ('N'): 32.62, + ('C','N'): -99.76, + ('Cl'): -24.31, + ('Br'): -35.42, + ('S'): -33.17, + ('Si'): -18.04 +} + +sigma_alpha = { + ('C'): 48.10, + ('O'): 193.06, + ('O','C'): 95.99, + ('C','N'): 152.20, + ('N','C'): -41.00 +} + +sigma_beta = { + ('C'): 32.87, + ('O'): -43.79, + ('O','O'):-128.16, + ('O','N'):79.13 +} + +# Molecular surface tension (cal/mol*AA^-2) +sigma_gamma = 0.35 +sigma_phi2 = -4.19 +sigma_psi2 = -6.68 +sigma_beta2 = 0.0 +gamma0 = 1.0 + +# rzz and delta r_zz in AA +r_zz = { + ('H','C'): [1.55, 0.3], + ('H','O'): [1.55, 0.3], + ('C','H'): [1.55, 0.3], + ('C','C'): [1.84, 0.3], + ('C','N'): [1.84, 0.3], + ('C','O'): [1.84, 0.3], + ('C','F'): [1.84, 0.3], + ('C','P'): [2.2, 0.3], + ('C','S'): [2.2, 0.3], + ('C','Cl'): [2.1, 0.3], + ('C','Br'): [2.3, 0.3], + ('C','I'): [2.6, 0.3], + ('N','C'): [1.84, 0.3], + ('N','C3'): [1.225, 0.065], + ('O','C'): [1.33, 0.1], + ('O','N'): [1.5, 0.3], + ('O','O'): [1.8, 0.3], + ('O','P'): [2.1, 0.3] +} + +def switch_function(R, r, dr): + return np.exp(dr/(R-dr-r)) if R