diff --git a/.gitattributes b/.gitattributes index 24ab84ebc..e6e9b020a 100644 --- a/.gitattributes +++ b/.gitattributes @@ -41,6 +41,10 @@ src/pedra/*.F90 licensefile=.githooks/LICENSE-Fortran src/pedra/*.f90 licensefile=.githooks/LICENSE-Fortran src/pedra/pedra_dlapack.f90 !licensefile +# tsless +src/tsless/*.F90 licensefile=.githooks/LICENSE-Fortran +src/tsless/*.f90 licensefile=.githooks/LICENSE-Fortran + # solver src/solver/*.hpp licensefile=.githooks/LICENSE-C++ src/solver/*.cpp licensefile=.githooks/LICENSE-C++ @@ -69,6 +73,7 @@ tests/iefpcm/*.cpp licensefile=.githooks/LICENSE-C++ tests/input/*.cpp licensefile=.githooks/LICENSE-C++ tests/numerical_quadrature/*.cpp licensefile=.githooks/LICENSE-C++ tests/utils/*.cpp licensefile=.githooks/LICENSE-C++ +tests/tsless/*.cpp licensefile=.githooks/LICENSE-C++ # tools doc/conf.py licensefile=.githooks/LICENSE-Python diff --git a/CHANGELOG.md b/CHANGELOG.md index d92ef66d3..bd111efd6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ ### Added +- A new cavity generator using C. S. Pomelli's [TsLess algorithm](http://dx.doi.org/10.1002/jcc.20076) - Green’s function for a spherical nanoparticle with **real** permittivity. The Green's function is known in analytical form from the work of Messina (J. Chem. Phys. 2002, 117 (24), 11062) and diff --git a/Pipfile b/Pipfile index c10830e68..7225c7559 100644 --- a/Pipfile +++ b/Pipfile @@ -7,14 +7,15 @@ name = "pypi" [packages] +PyYAML = "*" Pygments = "*" +Sphinx = "*" breathe = "*" docopt = "*" +fprettify = "*" matplotlib = "*" pyparsing = "*" -PyYAML = "*" recommonmark = "*" -Sphinx = "*" sphinx_rtd_theme = "*" sphinxcontrib-bibtex = "*" yapf = "*" diff --git a/Pipfile.lock b/Pipfile.lock index 32b70338d..d292f8af4 100644 --- a/Pipfile.lock +++ b/Pipfile.lock @@ -1,7 +1,7 @@ { "_meta": { "hash": { - "sha256": "48084d78a10e199dcab5f5012d4458b6c47ab8054af862ed5a9a900f79406c0d" + "sha256": "1daa1ebc9287a612c160b4189bf509c20f0488f69f6f6eb5a83caa8fc6b79162" }, "host-environment-markers": { "implementation_name": "cpython", @@ -9,9 +9,9 @@ "os_name": "posix", "platform_machine": "x86_64", "platform_python_implementation": "CPython", - "platform_release": "4.9.81", + "platform_release": "4.9.86", "platform_system": "Linux", - "platform_version": "#1-NixOS SMP Tue Feb 13 11:36:03 UTC 2018", + "platform_version": "#1-NixOS SMP Sat Mar 3 09:23:29 UTC 2018", "python_full_version": "3.6.4", "python_version": "3.6", "sys_platform": "linux" @@ -89,6 +89,20 @@ ], "version": "==0.14" }, + "fprettify": { + "hashes": [ + "sha256:438010a346c1dacad89ffac34ca8e0bb6c57ec07f8a7e8c01afc0f2f28843b2e", + "sha256:bdf565056aa35d134126c78a7af6ebbea7aa0658ab3b8ca623e5278537296f85", + "sha256:8a6452fad3ee4520cefe936cfadb762fe9e53aabead7896889eac4fedf173cc3" + ], + "version": "==0.3.1" + }, + "future": { + "hashes": [ + "sha256:e39ced1ab767b5936646cedba8bcce582398233d6a627067d4c6a454c90cfedb" + ], + "version": "==0.16.0" + }, "idna": { "hashes": [ "sha256:8c7309c718f94b3a625cb648ace320157ad16ff131ae0af362c9f21b80ef6ec4", @@ -110,6 +124,34 @@ ], "version": "==2.10" }, + "kiwisolver": { + "hashes": [ + "sha256:8b6a7b596ce1d2a6d93c3562f1178ebd3b7bb445b3b0dd33b09f9255e312a965", + "sha256:1a078f5dd7e99317098f0e0d490257fd0349d79363e8c923d5bb76428f318421", + "sha256:e0f910f84b35c36a3513b96d816e6442ae138862257ae18a0019d2fc67b041dc", + "sha256:aaec1cfd94f4f3e9a25e144d5b0ed1eb8a9596ec36d7318a504d813412563a85", + "sha256:f923406e6b32c86309261b8195e24e18b6a8801df0cfc7814ac44017bfcb3939", + "sha256:4329008a167fac233e398e8a600d1b91539dc33c5a3eadee84c0d4b04d4494fa", + "sha256:3b791ddf2aefc56382aadc26ea5b352e86a2921e4e85c31c1f770f527eb06ce4", + "sha256:379d97783ba8d2934d52221c833407f20ca287b36d949b4bba6c75274bcf6363", + "sha256:1aa0b55a0eb1bd3fa82e704f44fb8f16e26702af1a073cc5030eea399e617b56", + "sha256:ea36e19ac0a483eea239320aef0bd40702404ff8c7e42179a2d9d36c5afcb55c", + "sha256:0f7f532f3c94e99545a29f4c3f05637f4d2713e7fd91b4dd8abfc18340b86cd5", + "sha256:2874060b91e131ceeff00574b7c2140749c9355817a4ed498e82a4ffa308ecbc", + "sha256:95a25d9f3449046ecbe9065be8f8380c03c56081bc5d41fe0fb964aaa30b2195", + "sha256:79e5fe3ccd5144ae80777e12973027bd2f4f5e3ae8eb286cabe787bed9780138", + "sha256:9576cb63897fbfa69df60f994082c3f4b8e6adb49cccb60efb2a80a208e6f996", + "sha256:0ee4ed8b3ae8f5f712b0aa9ebd2858b5b232f1b9a96b0943dceb34df2a223bc3", + "sha256:66f82819ff47fa67a11540da96966fb9245504b7f496034f534b81cacf333861", + "sha256:b1c240d565e977d80c0083404c01e4d59c5772c977fae2c483f100567f50847b", + "sha256:53a5b27e6b5717bdc0125338a822605084054c80f382051fb945d2c0e6899a20", + "sha256:acb673eecbae089ea3be3dcf75bfe45fc8d4dcdc951e27d8691887963cf421c7", + "sha256:b15bc8d2c2848a4a7c04f76c9b3dc3561e95d4dabc6b4f24bfabe5fd81a0b14f", + "sha256:45813e0873bbb679334a161b28cb9606d9665e70561fd6caa8863e279b5e464b", + "sha256:ce3be5d520b4d2c3e5eeb4cd2ef62b9b9ab8ac6b6fedbaa0e39cdb6f50644278" + ], + "version": "==1.0.1" + }, "latexcodec": { "hashes": [ "sha256:3abc5cdbe4b6344c04a62db98bf30a2f3bc8f5bb9040296aafa39c8eeb737523", @@ -125,54 +167,54 @@ }, "matplotlib": { "hashes": [ - "sha256:31662334a4485455167072f80c57d2d2ef9ada4b93197b4851ceacee7269cb17", - "sha256:0c4a0cb10972b18049acde76e4611bcda0fa288a199a5f2c3c6f21ab208010cf", - "sha256:09ac1bde82673bfb4f1cb9ebd7fe258b29400f1be8914cb0891eca62a99cf7e6", - "sha256:295d2b135cb4d546fbb4b4d0b8d62d2132436adcf86221ece1d0bd4a3522128c", - "sha256:710bfe01633a4b734742163187bae99399de9f2c46cb6b40bd79b51467f7ba36", - "sha256:ae091625121ff5bf31515332ee604c90b55f9e6f6a02ac7160e5199f0422a211", - "sha256:fe8d2e29753fbbc1923e18e80487f09fdc01f270ff44da7156f147a3075876c0", - "sha256:9a87c1c8f195f0908b599268332608e5024b34cdd07375c68c15ac24b2a93a4e", - "sha256:ad987a871682d34d44a941cc47bf12f60e7866be7133498c432fb3d6fa3be651", - "sha256:59736af4bed61af336da60b8e97b700e695ffae9af08941c83e37ab0de3f151a", - "sha256:932ca62fb7c4edc377d6f0078995a1c004b7e4e9982a025e7022eb5eb1fc33ba", - "sha256:fb96735f76d5975e09eedc3502e9ff63dbe0c04aef312d0d6048a97f7993e667", - "sha256:2cc12ea94d80990e454a67ce829dc084bbdbfd0c95ed1f527dbd64688f184acd", - "sha256:97b9b56bbfce43deefbc2a089fa4afa8887e84041f17b090ac3d9f90a9e8e745", - "sha256:2452beef35840bdd9c3f613fc112ce6eb449cabca2bcb8a3f0d7eedd8e11d85d", - "sha256:adde3c9eb2145e5597a593877aea456e0bfe37f46cb3534caacce61b603c451a", - "sha256:fe5f8487aa872164a80b0c491f1110d4e2125391933caf17cae0a86df7950a72", - "sha256:725a3f12739d133adfa381e1b33bd70c6f64db453bfc536e148824816e568894" - ], - "version": "==2.1.2" + "sha256:3fd90b407d1ab0dae686a4200030ce305526ff20b85a443dc490d194114b2dfa", + "sha256:7b3d03c876684618e2a2be6abeb8d3a033c3a1bb38a786f751199753ef6227e6", + "sha256:abfd3d9390eb4f2d82cbcaa3a5c2834c581329b64eccb7a071ed9d5df27424f7", + "sha256:07055eb872fa109bd88f599bdb52065704b2e22d475b67675f345d75d32038a0", + "sha256:0f2f253d6d51f5ed52a819921f8a0a8e054ce0daefcfbc2557e1c433f14dc77d", + "sha256:1ef9fd285334bd6b0495b6de9d56a39dc95081577f27bafabcf28e0d318bed31", + "sha256:70f0e407fbe9e97f16597223269c849597047421af5eb8b60dbaca0382037e78", + "sha256:3fb2db66ef98246bafc04b4ef4e9b0e73c6369f38a29716844e939d197df816a", + "sha256:dc0ba2080fd0cfdd07b3458ee4324d35806733feb2b080838d7094731d3f73d9", + "sha256:4bb10087e09629ba3f9b25b6c734fd3f99542f93d71c5b9c023f28cb377b43a9", + "sha256:bc4d7481f0e8ec94cb1afc4a59905d6274b3b4c389aba7a2539e071766671735", + "sha256:4f6a516d5ef39128bb16af7457e73dde25c30625c4916d8fbd1cc7c14c55e691", + "sha256:8ff08eaa25c66383fe3b6c7eb288da3c22dcedc4b110a0b592b35f68d0e093b2", + "sha256:f26fba7fc68994ab2805d77e0695417f9377a00d36ba4248b5d0f1e5adb08d24", + "sha256:8944d311ce37bee1ba0e41a9b58dcf330ffe0cf29d7654c3d07c572215da68ac", + "sha256:45dac8589ef1721d7f2ab0f48f986694494dfcc5d13a3e43a5cb6c816276094e", + "sha256:9d12378d6a236aa38326e27f3a29427b63edce4ce325745785aec1a7535b1f85", + "sha256:4dc7ef528aad21f22be85e95725234c5178c0f938e2228ca76640e5e84d8cde8" + ], + "version": "==2.2.2" }, "numpy": { "hashes": [ - "sha256:e2335d56d2fd9fc4e3a3f2d3148aafec4962682375f429f05c45a64dacf19436", - "sha256:9b762e78739b6e021124adbea07611682db99cd3fca7f3c3a8b98b8f74ea5699", - "sha256:7d4c549e41507db4f04ec7cfab5597de8acf7871b16c9cf64cebcb9d39031ca6", - "sha256:b803306c4c201e7dcda0ce1b9a9c87f61a7c7ce43de2c60c8e56147b76849a1a", - "sha256:2da8dff91d489fea3e20155d41f4cd680de7d01d9a89fdd0ebb1bee6e72d3800", - "sha256:6b8c2daacbbffc83b4a2ba83a61aa3ce60c66340b07b962bd27b6c6bb175bee1", - "sha256:89b9419019c47ec87cf4cfca77d85da4611cc0be636ec87b5290346490b98450", - "sha256:49880b47d7272f902946dd995f346842c95fe275e2deb3082ef0495f0c718a69", - "sha256:3d7ddd5bdfb12ec9668edf1aa49a4a3eddb0db4661b57ea431477eb9a2468894", - "sha256:788e1757f8e409cd805a7cd82993cd9252fa19e334758a4c6eb5a8b334abb084", - "sha256:377def0873bbb1fbdedb14b3275b10a29b1b55619a3f7f775c4e7f9ce2461b9c", - "sha256:9501c9ccd081977ca5579a3ec4009d6baff6bacb04bf07214aade3324734195a", - "sha256:a1f5173df8190ef9c6235d260d70ca70c6fb029683ceb66e244c5cc6e335947a", - "sha256:12cf4b27039b88e407ad66894d99a957ef60fea0eeb442026af325add2ab264d", - "sha256:4e2fc841c8c642f7fd44591ef856ca409cedba6aea27928df34004c533839eee", - "sha256:e5ade7a69dccbd99c4fdbb95b6d091d941e62ffa588b0ed8fb0a2854118fef3f", - "sha256:6b1011ffc87d7e2b1b7bcc6dc21bdf177163658746ef778dcd21bf0516b9126c", - "sha256:a8bc80f69570e11967763636db9b24c1e3e3689881d10ae793cec74cf7a627b6", - "sha256:81b9d8f6450e752bd82e7d9618fa053df8db1725747880e76fb09710b57f78d0", - "sha256:e8522cad377cc2ef20fe13aae742cc265172910c98e8a0d6014b1a8d564019e2", - "sha256:a3d5dd437112292c707e54f47141be2f1100221242f07eda7bd8477f3ddc2252", - "sha256:c8000a6cbc5140629be8c038c9c9cdb3a1c85ff90bd4180ec99f0f0c73050b5e", - "sha256:fa0944650d5d3fb95869eaacd8eedbd2d83610c85e271bd9d3495ffa9bc4dc9c" - ], - "version": "==1.14.1" + "sha256:719d914f564f35cce4dc103808f8297c807c9f0297ac183ed81ae8b5650e698e", + "sha256:0f6a5ed0cd7ab1da11f5c07a8ecada73fc55a70ef7bb6311a4109891341d7277", + "sha256:d0928076d9bd8a98de44e79b1abe50c1456e7abbb40af7ef58092086f1a6c729", + "sha256:d858423f5ed444d494b15c4cc90a206e1b8c31354c781ac7584da0d21c09c1c3", + "sha256:20cac3123d791e4bf8482a580d98d6b5969ba348b9d5364df791ba3a666b660d", + "sha256:528ce59ded2008f9e8543e0146acb3a98a9890da00adf8904b1e18c82099418b", + "sha256:56e392b7c738bd70e6f46cf48c8194d3d1dd4c5a59fae4b30c58bb6ef86e5233", + "sha256:99051e03b445117b26028623f1a487112ddf61a09a27e2d25e6bc07d37d94f25", + "sha256:768e777cc1ffdbf97c507f65975c8686ebafe0f3dc8925d02ac117acc4669ce9", + "sha256:675e0f23967ce71067d12b6944add505d5f0a251f819cfb44bdf8ee7072c090d", + "sha256:a958bf9d4834c72dee4f91a0476e7837b8a2966dc6fcfc42c421405f98d0da51", + "sha256:bb370120de6d26004358611441e07acda26840e41dfedc259d7f8cc613f96495", + "sha256:f2b1378b63bdb581d5d7af2ec0373c8d40d651941d283a2afd7fc71184b3f570", + "sha256:a1413d06abfa942ca0553bf3bccaff5fdb36d55b84f2248e36228db871147dab", + "sha256:7f76d406c6b998d6410198dcb82688dcdaec7d846aa87e263ccf52efdcfeba30", + "sha256:a7157c9ac6bddd2908c35ef099e4b643bc0e0ebb4d653deb54891d29258dd329", + "sha256:0fd65cbbfdbf76bbf80c445d923b3accefea0fe2c2082049e0ce947c81fe1d3f", + "sha256:8c18ee4dddd5c6a811930c0a7c7947bf16387da3b394725f6063f1366311187d", + "sha256:0739146eaf4985962f07c62f7133aca89f3a600faac891ce6c7f3a1e2afe5272", + "sha256:07e21f14490324cc1160db101e9b6c1233c33985af4cb1d301dd02650fea1d7f", + "sha256:e6120d63b50e2248219f53302af7ec6fa2a42ed1f37e9cda2c76dbaca65036a7", + "sha256:6be6b0ca705321c178c9858e5ad5611af664bbdfae1df1541f938a840a103888", + "sha256:facc6f925c3099ac01a1f03758100772560a0b020fb9d70f210404be08006bcb" + ], + "version": "==1.14.2" }, "oset": { "hashes": [ @@ -182,10 +224,10 @@ }, "packaging": { "hashes": [ - "sha256:99276dc6e3a7851f32027a68f1095cd3f77c148091b092ea867a351811cfe388", - "sha256:5d50835fdf0a7edf0b55e311b7c887786504efea1177abd7e69329a8e5ea619e" + "sha256:e9215d2d2535d3ae866c3d6efc77d5b24a0192cce0ff20e42896cc0664f889c0", + "sha256:f019b770dd64e585a99714f1fd5e01c7a8f11b45635aa953fd41c689a657375b" ], - "version": "==16.8" + "version": "==17.1" }, "pybtex": { "hashes": [ @@ -221,10 +263,10 @@ }, "python-dateutil": { "hashes": [ - "sha256:95511bae634d69bc7329ba55e646499a842bc4ec342ad54a8cdb65645a0aad3c", - "sha256:891c38b2a02f5bb1be3e4793866c8df49c7d19baabf9c1bad62547e0b4866aca" + "sha256:3220490fb9741e2342e1cf29a503394fdac874bc39568288717ee67047ff29df", + "sha256:9d8074be4c993fbe4947878ce593052f71dac82932a677d49194d8ce9778002e" ], - "version": "==2.6.1" + "version": "==2.7.2" }, "pytz": { "hashes": [ @@ -289,10 +331,10 @@ }, "sphinx": { "hashes": [ - "sha256:41ae26acc6130ccf6ed47e5cca73742b80d55a134f0ab897c479bba8d3640b8e", - "sha256:da987de5fcca21a4acc7f67a86a363039e67ac3e8827161e61b91deb131c0ee8" + "sha256:7a606d77618a753adb79e13605166e3cf6a0e5678526e044236fc1ac43650910", + "sha256:5a1c9a0fec678c24b9a2f5afba240c04668edb7f45c67ce2ed008996b3f21ae2" ], - "version": "==1.7.1" + "version": "==1.7.2" }, "sphinx-rtd-theme": { "hashes": [ @@ -324,10 +366,10 @@ }, "yapf": { "hashes": [ - "sha256:0ba4765012218b67a8c67da5c045d1fb9171d786ffe8dce6c8d1e11cbbaba8eb", - "sha256:7f5efdb7edf0318b91e53721d934580a77153e24a222f52f6e1c3b7629aead43" + "sha256:dd23b52edbb4c0461d0383050f7886175b0df9ab8fd0b67edd41f94e25770993", + "sha256:7d8ae3567f3fb2d288f127d35e4decb3348c96cd091001e02e818465da618f90" ], - "version": "==0.20.2" + "version": "==0.21.0" } }, "develop": {} diff --git a/default.nix b/default.nix index f8b885b3e..08e37685e 100644 --- a/default.nix +++ b/default.nix @@ -3,8 +3,8 @@ let nixpkgs = (hostPkgs.fetchFromGitHub { owner = "NixOS"; repo = "nixpkgs-channels"; - rev = "nixos-unstable"; - sha256 = "0im8l87nghsp4z7swidgg2qpx9mxidnc0zs7a86qvw8jh7b6sbv2"; + rev = "nixos-18.03"; + sha256 = "1qiihxa2qdrij735rglfrx24nhn72kby44lkm1dxph395qy8fg1h"; }); in with import nixpkgs { diff --git a/requirements.txt b/requirements.txt index a3741f58d..4c40e448d 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,7 @@ Pygments breathe docopt +fprettify matplotlib pyparsing pyyaml diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 58a909841..f6b4c78b5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,26 +8,28 @@ include_directories( SYSTEM ${CMAKE_CURRENT_LIST_DIR}/utils/getkw ) +add_subdirectory(bi_operators) add_subdirectory(cavity) add_subdirectory(green) add_subdirectory(interface) add_subdirectory(metal) -add_subdirectory(bi_operators) add_subdirectory(pedra) add_subdirectory(solver) +add_subdirectory(tsless) add_subdirectory(utils) list(APPEND _objects - $ $<$:$> + $ + $ + $ $ $ $ - $ $ $ + $ $ - $ ) if(NOT STATIC_LIBRARY_ONLY) diff --git a/src/cavity/CMakeLists.txt b/src/cavity/CMakeLists.txt index 55af5f7bf..94236b931 100644 --- a/src/cavity/CMakeLists.txt +++ b/src/cavity/CMakeLists.txt @@ -6,6 +6,7 @@ Element.hpp GePolCavity.hpp ICavity.hpp RestartCavity.hpp +TsLessCavity.hpp ) # List of sources @@ -14,6 +15,7 @@ Element.cpp GePolCavity.cpp ICavity.cpp RestartCavity.cpp +TsLessCavity.cpp ) add_library(cavity OBJECT ${sources_list} ${headers_list}) diff --git a/src/cavity/Cavity.hpp b/src/cavity/Cavity.hpp index 78b011ef6..bcb261b1e 100644 --- a/src/cavity/Cavity.hpp +++ b/src/cavity/Cavity.hpp @@ -28,6 +28,7 @@ #include "GePolCavity.hpp" #include "ICavity.hpp" #include "RestartCavity.hpp" +#include "TsLessCavity.hpp" #include "utils/Factory.hpp" /*! @@ -51,6 +52,7 @@ inline Factory bootstrapFactory() { factory_.subscribe("GEPOL", createGePolCavity); factory_.subscribe("RESTART", createRestartCavity); + factory_.subscribe("TSLESS", createTsLessCavity); return factory_; } diff --git a/src/cavity/CavityData.hpp b/src/cavity/CavityData.hpp index 8328225c1..6d118dff1 100644 --- a/src/cavity/CavityData.hpp +++ b/src/cavity/CavityData.hpp @@ -42,6 +42,12 @@ struct CavityData { double area; /*! The radius of the spherical probe representing the solvent */ double probeRadius; + /*! The minimal distance between two sampling points on different spheres. Relevant + * for TsLessCavity */ + double minDistance; + /*! The maximum derivative order to be used in the definition of the smoothing + * function. Relevant for TsLessCavity */ + int derOrder; /*! Triggers the addition of spheres not centered on atoms, relevant for * GePolCavity */ double minimalRadius; @@ -52,12 +58,16 @@ struct CavityData { const Molecule & _molec, double _area, double _probeRadius, + double _minDistance, + int _derOrder, double _minRadius, const std::string & _fname) : cavityType(type), molecule(_molec), area(_area), probeRadius(_probeRadius), + minDistance(_minDistance), + derOrder(_derOrder), minimalRadius(_minRadius), filename(_fname) {} }; diff --git a/src/cavity/GePolCavity.cpp b/src/cavity/GePolCavity.cpp index f05a69393..b71cfce4d 100644 --- a/src/cavity/GePolCavity.cpp +++ b/src/cavity/GePolCavity.cpp @@ -83,7 +83,7 @@ void GePolCavity::build(const std::string & suffix, int maxsph, int maxvert) { - // This is a wrapper for the generatecavity_cpp function defined in the Fortran + // This is a wrapper for the pedra_driver function defined in the Fortran // code PEDRA. // Here we allocate the necessary arrays to be passed to PEDRA, in particular we // allow @@ -175,39 +175,39 @@ void GePolCavity::build(const std::string & suffix, int len_f_pedra = std::strlen(pedra.str().c_str()); // Go PEDRA, Go! TIMER_ON("GePolCavity::generatecavity_cpp"); - generatecavity_cpp(&maxts, - &maxsph, - &maxvert, - xtscor, - ytscor, - ztscor, - ar, - xsphcor, - ysphcor, - zsphcor, - rsph, - &nts, - &ntsirr, - &nSpheres_, - &addedSpheres, - xe, - ye, - ze, - rin, - mass, - &averageArea, - &probeRadius, - &minimalRadius, - &nr_gen, - &gen1, - &gen2, - &gen3, - nvert, - vert, - centr, - isphe, - pedra.str().c_str(), - &len_f_pedra); + pedra_driver(&maxts, + &maxsph, + &maxvert, + xtscor, + ytscor, + ztscor, + ar, + xsphcor, + ysphcor, + zsphcor, + rsph, + &nts, + &ntsirr, + &nSpheres_, + &addedSpheres, + xe, + ye, + ze, + rin, + mass, + &averageArea, + &probeRadius, + &minimalRadius, + &nr_gen, + &gen1, + &gen2, + &gen3, + nvert, + vert, + centr, + isphe, + pedra.str().c_str(), + &len_f_pedra); TIMER_OFF("GePolCavity::generatecavity_cpp"); // The "intensive" part of updating the spheres related class data members will be diff --git a/src/cavity/GePolCavity.hpp b/src/cavity/GePolCavity.hpp index 34cecef1e..442a7b067 100644 --- a/src/cavity/GePolCavity.hpp +++ b/src/cavity/GePolCavity.hpp @@ -94,18 +94,7 @@ class GePolCavity __final : public ICavity { void writeOFF(const std::string & suffix); }; -/*! \fn extern "C" void generatecavity_cpp(int * maxts, int * maxsph, int * maxvert, - * double * xtscor, double * ytscor, double * ztscor, - * double * ar, - * double * xsphcor, double * ysphcor, double * - * zsphcor, double * rsph, - * int * nts, int * ntsirr, int * nesfp, int * - * addsph, - * double * xe, double * ye, double * ze, double * - * rin, - * double * avgArea, double * rsolv, double * ret, - * int * nr_gen, int * gen1, int * gen2, int * gen3, - * int * nvert, double * vert, double * centr) +/*! \brief Interface to the Fortran PEDRA code * \param[in] maxts maximum number of tesserae allowed * \param[in] maxsph maximum number of spheres allowed * \param[in] maxvert maximum number of vertices allowed @@ -144,39 +133,39 @@ class GePolCavity __final : public ICavity { * \param[out] vert coordinates of tesserae vertices * \param[out] centr centers of arcs defining the edges of the tesserae */ -extern "C" void generatecavity_cpp(int * maxts, - int * maxsph, - int * maxvert, - double * xtscor, - double * ytscor, - double * ztscor, - double * ar, - double * xsphcor, - double * ysphcor, - double * zsphcor, - double * rsph, - int * nts, - int * ntsirr, - int * nesfp, - int * addsph, - double * xe, - double * ye, - double * ze, - double * rin, - double * masses, - double * avgArea, - double * rsolv, - double * ret, - int * nr_gen, - int * gen1, - int * gen2, - int * gen3, - int * nvert, - double * vert, - double * centr, - int * isphe, - const char * pedra, - int * len_f_pedra); +extern "C" void pedra_driver(int * maxts, + int * maxsph, + int * maxvert, + double * xtscor, + double * ytscor, + double * ztscor, + double * ar, + double * xsphcor, + double * ysphcor, + double * zsphcor, + double * rsph, + int * nts, + int * ntsirr, + int * nesfp, + int * addsph, + double * xe, + double * ye, + double * ze, + double * rin, + double * masses, + double * avgArea, + double * rsolv, + double * ret, + int * nr_gen, + int * gen1, + int * gen2, + int * gen3, + int * nvert, + double * vert, + double * centr, + int * isphe, + const char * pedra, + int * len_f_pedra); ICavity * createGePolCavity(const CavityData & data); } // namespace cavity diff --git a/src/cavity/TsLessCavity.cpp b/src/cavity/TsLessCavity.cpp new file mode 100644 index 000000000..9f667c4ee --- /dev/null +++ b/src/cavity/TsLessCavity.cpp @@ -0,0 +1,341 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PCMSolver is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "TsLessCavity.hpp" + +#include +#include +#include +#include +#include + +#include "Config.hpp" + +#include + +#include "CavityData.hpp" +#include "utils/Sphere.hpp" +#include "utils/Symmetry.hpp" + +namespace pcm { +namespace cavity { +TsLessCavity::TsLessCavity(const Molecule & molec, + double a, + double pr, + double minR, + double minD, + int der) + : ICavity(molec), + averageArea_(a), + probeRadius_(pr), + minimalRadius_(minR), + minDistance_(minD), + derOrder_(der) { + TIMER_ON("TsLessCavity::build from Molecule"); + build(10000, 200, 25000); + TIMER_OFF("TsLessCavity::build from Molecule"); +} + +TsLessCavity::TsLessCavity(const Sphere & sph, + double a, + double pr, + double minR, + double minD, + int der) + : ICavity(sph), + averageArea_(a), + probeRadius_(pr), + minimalRadius_(minR), + minDistance_(minD), + derOrder_(der) { + TIMER_ON("TsLessCavity::build from a single sphere"); + build(10000, 200, 25000); + TIMER_OFF("TsLessCavity::build from a single sphere"); +} + +TsLessCavity::TsLessCavity(const std::vector & sph, + double a, + double pr, + double minR, + double minD, + int der) + : ICavity(sph), + averageArea_(a), + probeRadius_(pr), + minimalRadius_(minR), + minDistance_(minD), + derOrder_(der) { + TIMER_ON("TsLessCavity::build from list of spheres"); + build(10000, 200, 25000); + TIMER_OFF("TsLessCavity::build from list of spheres"); +} + +void TsLessCavity::build(int maxts, int maxsph, int maxvert) { + // This is a wrapper for the generatecavity_cpp_ function defined in the Fortran + // code TsLess. Here we allocate the necessary arrays to be passed to TsLess, in + // particular we allow for the insertion of additional spheres as in the most + // general formulation of the GePol algorithm. + + int lwork = maxts * maxsph; + double * xtscor = new double[maxts]; + double * ytscor = new double[maxts]; + double * ztscor = new double[maxts]; + double * ar = new double[maxts]; + double * xsphcor = new double[maxts]; + double * ysphcor = new double[maxts]; + double * zsphcor = new double[maxts]; + double * rsph = new double[maxts]; + double * work = new double[lwork]; + + // Clean-up possible heap-crap + std::fill_n(xtscor, maxts, 0.0); + std::fill_n(ytscor, maxts, 0.0); + std::fill_n(ztscor, maxts, 0.0); + std::fill_n(ar, maxts, 0.0); + std::fill_n(xsphcor, maxts, 0.0); + std::fill_n(ysphcor, maxts, 0.0); + std::fill_n(zsphcor, maxts, 0.0); + std::fill_n(rsph, maxts, 0.0); + std::fill_n(work, lwork, 0.0); + + int nts = 0; + int ntsirr = 0; + + // If there's an overflow in the number of spheres TsLess will die. + // The maximum number of spheres in PEDRA is set to 200 (primitive+additional) + // so the integer here declared is just to have enough space C++ side to pass + // everything back. + int maxAddedSpheres = 200; + + // maximum number of spheres we allow the algorithm to add to our original set. + // If this number is exceeded, then the algorithm crashes (should look into + // this...) After the cavity is generated we will update ALL the class data + // members, both related to spheres and finite elements so that the cavity is fully + // formed. + + Eigen::VectorXd xv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres); + Eigen::VectorXd yv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres); + Eigen::VectorXd zv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres); + Eigen::VectorXd radii_scratch = + Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres); // Not to be confused with + // the data member + // inherited from Cavity!!! + + for (int i = 0; i < nSpheres_; ++i) { + for (int j = 0; j < 3; ++j) { + xv(i) = sphereCenter_(0, i); + yv(i) = sphereCenter_(1, i); + zv(i) = sphereCenter_(2, i); + } + radii_scratch(i) = sphereRadius_(i); + } + + double * xe = xv.data(); + double * ye = yv.data(); + double * ze = zv.data(); + + double * rin = radii_scratch.data(); + double * mass = new double[molecule_.nAtoms()]; + for (size_t i = 0; i < molecule_.nAtoms(); ++i) { + mass[i] = molecule_.masses(i); + } + + addedSpheres = 0; + // Number of generators and generators of the point group + int nr_gen = molecule_.pointGroup().nrGenerators(); + int gen1 = molecule_.pointGroup().generators(0); + int gen2 = molecule_.pointGroup().generators(1); + int gen3 = molecule_.pointGroup().generators(2); + + int weightFunction = 1; + // Go TsLess, Go! + TIMER_ON("TsLessCavity::tsless_driver"); + tsless_driver(&maxts, + &maxsph, + &maxvert, + &nSpheres_, + &nts, + &ntsirr, + &addedSpheres, + xtscor, + ytscor, + ztscor, + ar, + xsphcor, + ysphcor, + zsphcor, + rsph, + xe, + ye, + ze, + rin, + mass, + &nr_gen, + &gen1, + &gen2, + &gen3, + &averageArea_, + &minDistance_, + &derOrder_, + &weightFunction, + &probeRadius_, + work); + TIMER_OFF("TsLessCavity::tsless_driver"); + + // The "intensive" part of updating the spheres related class data members will be + // of course executed iff addedSpheres != 0 + if (addedSpheres != 0) { + // Save the number of original spheres + int orig = nSpheres_; + // Update the nSpheres + nSpheres_ += addedSpheres; + // Resize sphereRadius and sphereCenter... + sphereRadius_.resize(nSpheres_); + sphereCenter_.resize(Eigen::NoChange, nSpheres_); + // Transfer radii and centers. + // Eigen has no push_back function, so we need to traverse all the spheres... + sphereRadius_ = radii_scratch.head(nSpheres_); + for (int i = 0; i < nSpheres_; ++i) { + sphereCenter_(0, i) = xv(i); + sphereCenter_(1, i) = yv(i); + sphereCenter_(2, i) = zv(i); + } + // Now grow the vector containing the list of spheres + for (int i = orig; i < nSpheres_; ++i) { + spheres_.push_back(Sphere(sphereCenter_.col(i), sphereRadius_(i))); + } + } + + nElements_ = static_cast(nts); + nIrrElements_ = static_cast(ntsirr); + elementCenter_.resize(Eigen::NoChange, nElements_); + elementSphereCenter_.resize(Eigen::NoChange, nElements_); + elementNormal_.resize(Eigen::NoChange, nElements_); + elementArea_.resize(nElements_); + elementRadius_.resize(nElements_); + for (int i = 0; i < nElements_; ++i) { + elementCenter_(0, i) = xtscor[i]; + elementCenter_(1, i) = ytscor[i]; + elementCenter_(2, i) = ztscor[i]; + elementArea_(i) = ar[i]; + elementSphereCenter_(0, i) = xsphcor[i]; + elementSphereCenter_(1, i) = ysphcor[i]; + elementSphereCenter_(2, i) = zsphcor[i]; + elementRadius_(i) = rsph[i]; + } + + elementNormal_ = elementCenter_ - elementSphereCenter_; + for (int i = 0; i < nElements_; ++i) { + elementNormal_.col(i) /= elementNormal_.col(i).norm(); + } + + // Fill elements_ vector + for (int i = 0; i < nElements_; ++i) { + bool irr = false; + int nv = 1; // TsLess does not generate spherical polygons!! + // TsLess puts the irreducible tesserae first (? Check with Cris!) + if (i < nIrrElements_) + irr = true; + Sphere sph(elementSphereCenter_.col(i), elementRadius_(i)); + Eigen::Matrix3Xd vertices, arcs; + vertices.resize(Eigen::NoChange, nv); + arcs.resize(Eigen::NoChange, nv); + // FIXME index of the sphere the element belongs to + elements_.push_back(Element(nv, + 0, + elementArea_(i), + elementCenter_.col(i), + elementNormal_.col(i), + irr, + sph, + vertices, + arcs)); + } + + delete[] xtscor; + delete[] ytscor; + delete[] ztscor; + delete[] ar; + delete[] xsphcor; + delete[] ysphcor; + delete[] zsphcor; + delete[] rsph; + delete[] work; + delete[] mass; + + built = true; +} + +std::ostream & TsLessCavity::printCavity(std::ostream & os) { + os << "Cavity type: TsLess" << std::endl; + os << "Average point weight = " << averageArea_ << " AU^2" << std::endl; + os << "Minimal distance between sampling points = " << minDistance_ << " AU" + << std::endl; + os << "Switch function is of class C^" << derOrder_ << std::endl; + os << "Addition of extra spheres enabled" << std::endl; + os << "Probe radius = " << probeRadius_ << " AU" << std::endl; + os << "Number of spheres = " << nSpheres_ + << " [initial = " << nSpheres_ - addedSpheres << "; added = " << addedSpheres + << "]" << std::endl; + os << "Number of finite elements = " << nElements_; + os << "Number of irreducible finite elements = " << nIrrElements_ << std::endl; + os << "============ Spheres list (in Angstrom)" << std::endl; + os << " Sphere on Radius Alpha X Y Z \n"; + os << "-------- ---- -------- ------- ----------- ----------- -----------\n"; + // Print original set of spheres + int original = nSpheres_ - addedSpheres; + Eigen::IOFormat CleanFmt(6, Eigen::DontAlignCols, " ", "\n", "", ""); + for (int i = 0; i < original; ++i) { + os << std::setw(4) << i + 1; + os << " " << molecule_.atoms()[i].symbol << " "; + os << std::fixed << std::setprecision(4) << molecule_.atoms()[i].radius; + os << std::fixed << std::setprecision(2) << " " + << molecule_.atoms()[i].radiusScaling << " "; + os << (molecule_.geometry().col(i).transpose() * bohrToAngstrom()) + .format(CleanFmt); + os << std::endl; + } + // Print added spheres + for (int j = 0; j < addedSpheres; ++j) { + int idx = original + j; + os << std::setw(4) << idx + 1; + os << " Du "; + os << std::fixed << std::setprecision(4) + << sphereRadius_(idx) * bohrToAngstrom(); + os << std::fixed << std::setprecision(2) << " 1.00"; + os << (sphereCenter_.col(idx).transpose() * bohrToAngstrom()).format(CleanFmt); + os << std::endl; + } + return os; +} + +ICavity * createTsLessCavity(const CavityData & data) { + return new TsLessCavity(data.molecule, + data.area, + data.probeRadius, + data.minimalRadius, + data.minDistance, + data.derOrder); +} +} // namespace cavity +} // namespace pcm diff --git a/src/cavity/TsLessCavity.hpp b/src/cavity/TsLessCavity.hpp new file mode 100644 index 000000000..5a65f9a5f --- /dev/null +++ b/src/cavity/TsLessCavity.hpp @@ -0,0 +1,158 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PCMSolver is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#pragma once + +#include +#include +#include + +#include "Config.hpp" + +namespace pcm { +struct CavityData; +} // namespace pcm + +#include "ICavity.hpp" + +/*! \file TsLessCavity.hpp */ + +namespace pcm { +namespace cavity { +/*! \class TsLessCavity + * \brief A class for TsLess cavity. + * \author Roberto Di Remigio + * \date 2013 + * + * This class is a wrapper for the Fortran routines generating + * a tessellationless grid on the cavity surface. The original + * algoritm is described in \cite Pomelli2004 + */ +class TsLessCavity __final : public ICavity { +public: + TsLessCavity() {} + TsLessCavity(const Molecule & molec, + double a, + double pr, + double minR, + double minD, + int der); + TsLessCavity(const Sphere & sph, + double a, + double pr, + double minR, + double minD, + int der); + TsLessCavity(const std::vector & sph, + double a, + double pr, + double minR, + double minD, + int der); + virtual ~TsLessCavity() {} + friend std::ostream & operator<<(std::ostream & os, TsLessCavity & cavity) { + return cavity.printCavity(os); + } + +private: + double averageArea_; + double probeRadius_; + double minimalRadius_; + double minDistance_; + int derOrder_; + int addedSpheres; + virtual std::ostream & printCavity(std::ostream & os) __override; + virtual void makeCavity() __override { build(10000, 200, 25000); } + /*! \brief Driver for TsLess Fortran module. + * \param[in] maxts maximum number of tesserae + * \param[in] maxsp maximum number of spheres (original + added) + * \param[in] maxvert maximum number of vertices + */ + void build(int maxts, int maxsp, int maxvert); +}; + +/*! \brief Fortran interface function to TsLess cavity generation + * \param[in] maxts maximum number of tesserae allowed + * \param[in] maxsph maximum number of spheres allowed + * \param[in] maxvert maximum number of vertices allowed + * \param[out] nesfp number of spheres (original + added) + * \param[out] nts number of generated tesserae + * \param[out] ntsirr number of generated irreducible tesserae + * \param[out] addsph number of added spheres + * \param[out] xtscor x-coordinate of tesserae centers (dimension maxts) + * \param[out] ytscor y-coordinate of tesserae centers (dimension maxts) + * \param[out] ztscor z-coordinate of tesserae centers (dimension maxts) + * \param[out] ar area of the tessera (dimension maxts) + * \param[out] xsphcor x-coordinate of the sphere center the tessera belongs to + * (dimension maxts) \param[out] ysphcor y-coordinate of the sphere center the + * tessera belongs to (dimension maxts) \param[out] zsphcor z-coordinate of the + * sphere center the tessera belongs to (dimension maxts) \param[out] rsph radii of + * the sphere the tessera belongs to, i.e. its curvature (dimension maxts) + * \param[out] xe x-coordinate of the sphere center (dimension nSpheres_ + + * maxAddedSpheres) \param[out] ye y-coordinate of the sphere center (dimension + * nSpheres_ + maxAddedSpheres) \param[out] ze z-coordinate of the sphere center + * (dimension nSpheres_ + maxAddedSpheres) \param[out] rin radius of the spheres + * (dimension nSpheres_ + maxAddedSpheres) \param[in] masses atomic masses (for + * inertia tensor formation in TSLESS) \param[in] nr_gen number of symmetry + * generators \param[in] gen1 first generator \param[in] gen2 second generator + * \param[in] gen3 third generator + * \param[in] avgArea average tesserae area + * \param[in] dmin mininal distance between sampling points + * \param[in] nord maximum order of continuous derivative of weight function + * \param[in] ifun whether to use the normalized or unnormalized form of the weight + * function \param[in] rsolv solvent probe radius \param[in] work scratch space + */ +extern "C" void tsless_driver(int * maxts, + int * maxsph, + int * maxvert, + int * nesfp, + int * nts, + int * ntsirr, + int * addsph, + double * xtscor, + double * ytscor, + double * ztscor, + double * ar, + double * xsphcor, + double * ysphcor, + double * zsphcor, + double * rsph, + double * xe, + double * ye, + double * ze, + double * rin, + double * masses, + int * nr_gen, + int * gen1, + int * gen2, + int * gen3, + double * avgArea, + double * dmin, + int * nord, + int * ifun, + double * rsolv, + double * work); + +ICavity * createTsLessCavity(const CavityData & data); +} // namespace cavity +} // namespace pcm diff --git a/src/interface/Input.cpp b/src/interface/Input.cpp index 98d7c5ce8..8ae46d968 100644 --- a/src/interface/Input.cpp +++ b/src/interface/Input.cpp @@ -73,6 +73,8 @@ void Input::reader(const std::string & filename) { cavityType_ = cavity.getStr("TYPE"); area_ = cavity.getDbl("AREA"); + minDistance_ = cavity.getDbl("MINDISTANCE"); + derOrder_ = cavity.getInt("DERORDER"); if (cavityType_ == "RESTART") { cavFilename_ = cavity.getStr("NPZFILE"); } @@ -201,6 +203,8 @@ void Input::reader(const PCMInput & host_input) { cavityType_ = detail::trim_and_upper(host_input.cavity_type); area_ = host_input.area * angstrom2ToBohr2(); + minDistance_ = host_input.min_distance * angstromToBohr(); + derOrder_ = host_input.der_order; if (cavityType_ == "RESTART") { cavFilename_ = detail::trim(host_input.restart_name); // No case conversion here! } @@ -339,8 +343,14 @@ void Input::initMolecule() { } CavityData Input::cavityParams() const { - return CavityData( - cavityType_, molecule_, area_, probeRadius_, minimalRadius_, cavFilename_); + return CavityData(cavityType_, + molecule_, + area_, + probeRadius_, + minDistance_, + derOrder_, + minimalRadius_, + cavFilename_); } GreenData Input::insideGreenParams() const { diff --git a/src/interface/Input.hpp b/src/interface/Input.hpp index e2d17ada8..9a7af66e0 100644 --- a/src/interface/Input.hpp +++ b/src/interface/Input.hpp @@ -170,8 +170,12 @@ class PCMSolver_EXPORT Input { std::string cavityType_; /// Filename for the .npz cavity restart file std::string cavFilename_; - /// GePol cavity average element area + /// GePol and TsLess cavity average element area double area_; + /// TsLess cavity minimal distance between sampling points + double minDistance_; + /// TsLess cavity maximum derivative order of switch function + int derOrder_; /// Whether the radii should be scaled by 1.2 bool scaling_; /// The set of radii to be used diff --git a/src/pedra/pedra_cavity_interface.F90 b/src/pedra/pedra_cavity_interface.F90 index 07d962396..e88eb723b 100644 --- a/src/pedra/pedra_cavity_interface.F90 +++ b/src/pedra/pedra_cavity_interface.F90 @@ -27,14 +27,14 @@ ! ! RDR, 280114. Put things in makecav.F inside here directly. ! -subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & +subroutine pedra_driver(maxts_, maxsph_, maxvert_, & xtscor_, ytscor_, ztscor_, ar_, & xsphcor_, ysphcor_, zsphcor_, rsph_, & nts_, ntsirr_, nesfp_, addsph_, & xe_, ye_, ze_, rin_, masses_, avgArea_, rsolv_, ret_, & nr_gen_, gen1_, gen2_, gen3_, & nvert_, vert_, centr_, isphe_, pedra_, len_pedra_) & - bind(C, name='generatecavity_cpp') + bind(C, name='pedra_driver') use, intrinsic :: iso_c_binding use pedra_precision @@ -70,7 +70,6 @@ subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & character(len=len_pedra_) :: pedra type(point_group) :: pgroup - pedra = carray_to_fstring(pedra_) !lvpri = 121201_regint_k ! The following INQUIRE statement returns whether the file named cavity.off is @@ -172,4 +171,4 @@ subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & close(pedra_unit) -end subroutine generatecavity_cpp +end subroutine diff --git a/src/tsless/CMakeLists.txt b/src/tsless/CMakeLists.txt new file mode 100644 index 000000000..9803c4f21 --- /dev/null +++ b/src/tsless/CMakeLists.txt @@ -0,0 +1,20 @@ +# List of sources +list(APPEND sources_list +tsless_cavity.f90 +tsless_lapack.f90 +tsless_precision.f90 +tsless_symmetry.f90 +tsless_weight_function.f90 +) + +add_library(tsless OBJECT ${sources_list}) +set_target_properties(tsless + PROPERTIES + POSITION_INDEPENDENT_CODE 1 + CXX_VISIBILITY_PRESET hidden + VISIBILITY_INLINES_HIDDEN 1 + ) +target_compile_definitions(tsless + PUBLIC + PCMSolver_EXPORTS + ) diff --git a/src/tsless/tsless_cavity.f90 b/src/tsless/tsless_cavity.f90 new file mode 100644 index 000000000..f773ca5ee --- /dev/null +++ b/src/tsless/tsless_cavity.f90 @@ -0,0 +1,390 @@ +! +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. +! +! This file is part of PCMSolver. +! +! PCMSolver is free software: you can redistribute it and/or modify +! it under the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! PCMSolver is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with PCMSolver. If not, see . +! +! For information on the complete list of contributors to the +! PCMSolver API, see: +! + +!> Originally written by Christian Silvio Pomelli (ca. 2003-2004) +!> Introduced in PCMSolver by Christian Silvio Pomelli (ca. 2014) +!> Wrapped in a F90 module by Roberto Di Remigio (2015) +module tsless_cavity + +use, intrinsic :: iso_c_binding +use tsless_precision +use tsless_symmetry +use tsless_weight_function + +implicit none + +public tsless_driver + +private + +contains + + !> \brief Interface function to TsLess cavity generation + !> \author Roberto Di Remigio, Christian S. Pomelli + !> \year 2014, 2015 + !> \param[in] maxts maximum number of tesserae allowed + !> \param[in] maxsph maximum number of spheres allowed + !> \param[in] maxvert maximum number of vertices allowed + !> \param[out] nesfp number of spheres (original + added) + !> \param[out] nts number of generated tesserae + !> \param[out] ntsirr number of generated irreducible tesserae + !> \param[out] addsph number of added spheres + !> \param[out] xtscor x-coordinate of tesserae centers + !> \param[out] ytscor y-coordinate of tesserae centers + !> \param[out] ztscor z-coordinate of tesserae centers + !> \param[out] ar area of the tessera + !> \param[out] xsphcor x-coordinate of the sphere center the tessera belongs to + !> \param[out] ysphcor y-coordinate of the sphere center the tessera belongs to + !> \param[out] zsphcor z-coordinate of the sphere center the tessera belongs to + !> \param[out] rsph radii of the sphere the tessera belongs to, i.e. its curvature + !> \param[out] xe x-coordinate of the sphere center + !> \param[out] ye y-coordinate of the sphere center + !> \param[out] ze z-coordinate of the sphere center + !> \param[out] rin radius of the spheres + !> \param[in] masses atomic masses (for inertia tensor formation in TSLESS) + !> \param[in] nr_gen number of symmetry generators + !> \param[in] gen1 first generator + !> \param[in] gen2 second generator + !> \param[in] gen3 third generator + !> \param[in] avgArea average tesserae area + !> \param[in] dmin mininal distance between sampling points + !> \param[in] nord maximum order of continuous derivative of weight function + !> \param[in] ifun whether to use the normalized or unnormalized form of the weight function + !> \param[in] rsolv solvent probe radius + !> \param[in] work scratch space + subroutine tsless_driver(maxts, maxsph, maxvert, nesfp, nts, ntsirr, addsph, & + xtscor, ytscor, ztscor, ar, & + xsphcor, ysphcor, zsphcor, rsph, & + xe, ye, ze, rin, & + masses, & + nr_gen, gen1, gen2, gen3, & + avgArea, dmin, nord, ifun, rsolv, work) & + bind(C, name='tsless_driver') + + !> Passed variables + integer(c_size_t), intent(in) :: maxts, maxsph, maxvert + integer(c_int), intent(in) :: nesfp + real(c_double), intent(out) :: xtscor(maxts), ytscor(maxts), ztscor(maxts), ar(maxts) + real(c_double), intent(out) :: xsphcor(maxts), ysphcor(maxts), zsphcor(maxts), rsph(maxts) + real(c_double), intent(in) :: xe(nesfp), ye(nesfp), ze(nesfp), rin(nesfp) + real(c_double), intent(in) :: masses(nesfp) + real(c_double), intent(in) :: work(maxts*maxsph) + real(c_double), intent(in) :: avgArea, rsolv, dmin + integer(c_int), intent(in) :: nord, ifun + integer(c_int), intent(out) :: nts, ntsirr, addsph + integer(c_int), intent(in) :: nr_gen, gen1, gen2, gen3 + !> Parameters + real(kind=dp), parameter :: pi = acos(-1.0_dp) + real(kind=dp), parameter :: tpi = 2.0_dp * pi + real(kind=dp), parameter :: fpi = 4.0_dp * pi + !> Local variables + integer :: print_unit + integer :: i, j, k, ntssp + logical :: tsless_file_exists, tsless_file_opened + integer(kind=regint_k) :: nesf_orig + type(weight_function) :: weights + real(kind=dp) :: x, y, z, r, dscale(3) + real(kind=dp) :: dist, dist2, rho, t, w0 + type(point_group) :: pgroup + + print_unit = 121221_regint_k + tsless_file_exists = .false. + inquire(file = 'TSLESS.OUT', exist = tsless_file_exists, opened = tsless_file_opened) + if (.not. tsless_file_opened) then + open(print_unit, & + file = 'TSLESS.OUT', & + status = 'replace', & + form = 'formatted', & + access = 'sequential') + end if + + !> Save original number of spheres + nesf_orig = nesfp + !> Create point group + pgroup = build_point_group(nr_gen, gen1, gen2, gen3, print_unit) + write(print_unit, *) 'Passed from host: nord, dmin, avgarea', nord, dmin, avgarea + + ! version 0.1 10/1/2014: plain implementation + ! added sphere positions, radii and derivatives (TBI) + ! generation and weighting of sampling points + + ! Generate weight function + weights = create_weight_function(nderiv=nord, normalization=ifun, dmin=dmin, printer=print_unit) + nts = 0_regint_k + + ! main loop on spheres (in future will implement linear scaling) + LoopOnSpheres: do i = 1, nesfp + ! point generation with weighting according to Leopardi scheme [ref] + ntssp = int(fpi * rin(i)**2 / avgArea) + 1_regint_k + !> Original weight factor + w0 = fpi * rin(i)**2 / ntssp + !> Interaction radius + rho = sqrt(w0 * pi) + write(print_unit, '(a, f12.5, a, f12.5)') 'Original weight factor ', w0, ' and interaction radius ', rho + x = xe(i) + y = ye(i) + z = ze(i) + r = rin(i) + + ! generate sphere sampling points + write(print_unit, '(a, i4, a, i4, a)') 'Generating ', ntssp, ' sampling points on sphere', i, ' according to Leopardi' + write(print_unit, '(a, f12.5, a, f12.5, a, f12.5, a, f12.5)') 'Center = (', x, ', ', y, ', ', z, ') Radius = ', r + call leopardi_grid(ntssp, nts, x, y, z, r, xtscor, ytscor, ztscor, print_unit) + write(print_unit, '(a)') 'Sampling points generated, now refining for intersections between spheres' + + ! loop on sphere sampling points + LoopOnSamplingPoints: do j = nts + 1, nts + ntssp + ar(j) = w0 + ! loop on intersecting spheres + LoopOnIntersectingSpheres: do k = 1, nesfp + if (k .eq. i) then + goto 10 + endif + dist2 = (xe(k) - xtscor(j))**2 + (ye(k) - ytscor(j))**2 + (ze(k) - ztscor(j))**2 + dist = sqrt(dist2) + t = (dist - rin(k)) / rho + write(print_unit, *) 't = ', t + if (t .le. dmin) then + ar(j) = -1.0_dp + write(print_unit, *) 'point ', j, ' was deleted' + goto 20 + else if(t .le. 1.0_dp) then + write(print_unit, *) 'point ', j, ' was scaled' + write(print_unit, *) 'Original value ', ar(j) + dscale = scaling_factor(weights, t, dmin) + ar(j) = ar(j) * dscale(1) + write(print_unit, *) 'Current value ', ar(j), ' Scaling factor ', dscale(1) + endif + 10 continue + end do LoopOnIntersectingSpheres + 20 continue + + ! save points with non zero weights + SavePoints: if (ar(j) .gt. 0.0_dp) then + nts = nts + 1_regint_k + xtscor(nts) = xtscor(j) + ytscor(nts) = ytscor(j) + ztscor(nts) = ztscor(j) + ar(nts) = ar(j) + xsphcor(nts) = x + ysphcor(nts) = y + zsphcor(nts) = z + rsph(nts) = r + end if SavePoints + end do LoopOnSamplingPoints + end do LoopOnSpheres + + !> Number of irreducible sampling points, hardcoded to nts for the moment + ntsirr = nts + !> Number of added spheres + addsph = nesfp - nesf_orig + + call destroy_weight_function(weights) + + write(print_unit, *) ' Surface = ', cavity_surface(nts, ar), ' volume = ', cavity_volume(nts, rsph, ar) + close(print_unit) + + call write_visualization_file(nts, nesfp, xtscor, ytscor, ztscor, ar, xe, ye, ze, rin) + + end subroutine tsless_driver + + !> \brief Create and translate Leopardi points + !> \author Christian S. Pomelli + !> \param[in] ntssp number of points on the sphere + !> \param[in] nts a counter for the total number of points generated + !> \param[in] x x-coordinate of sphere center + !> \param[in] y y-coordinate of sphere center + !> \param[in] z z-coordinate of sphere center + !> \param[in] r radius of sphere + !> \param[in] xtscor x-coordinates of points on the sphere + !> \param[in] ytscor y-coordinates of points on the sphere + !> \param[in] ztscor z-coordinates of points on the sphere + !> \param[in] print_unit + subroutine leopardi_grid(ntssp, nts, x, y, z, r, xtscor, ytscor, ztscor, print_unit) + + use, intrinsic :: iso_fortran_env, only: output_unit + + !> Passed variables + integer(kind=regint_k), intent(in) :: ntssp, nts + integer, optional, intent(in) :: print_unit + real(kind=dp), intent(in) :: x, y, z, r + real(kind=dp), intent(inout) :: xtscor(*), ytscor(*), ztscor(*) + !> Parameters + real(kind=dp), parameter :: pi = acos(-1.0_dp) + real(kind=dp), parameter :: tpi = 2.0_dp * pi + real(kind=dp), parameter :: fpi = 4.0_dp * pi + !> Local variables + integer :: i, j, k, n, print_out + real(kind=dp) :: thetaf(ntssp), ypsi(ntssp), theta(ntssp), apsi(ntssp), m(ntssp) + real(kind=dp) :: vr, thetac, delta1, deltaf + real(kind=dp) :: thetam, costm, sintm, omegaj + + if (present(print_unit)) then + print_out = print_unit + else + print_out = output_unit + end if + + thetaf = 0.0_dp + ypsi = 0.0_dp + theta = 0.0_dp + apsi = 0.0_dp + m = 0.0_dp + ! Notation consistent with [ref] + vr = fpi / ntssp + thetac = 2.0_dp * (asin(sqrt(vr / fpi))) + delta1 = sqrt(vr) + + ! Polar points + ! North pole + xtscor(nts+1) = x + ytscor(nts+1) = y + ztscor(nts+1) = z + r + ! South pole + xtscor(nts+ntssp) = x + ytscor(nts+ntssp) = y + ztscor(nts+ntssp) = z - r + + ! two points case (rare) + if (ntssp == 2) return + + ! Number of collars + n = int((pi - 2.0_dp * thetac) / delta1 + 0.5_dp) + if (n == 0) n = 1 + deltaf = (pi - 2.0_dp * thetac) / (1.0_dp * n) + + ! Colatitudes of collars + CollarColatitudes: do i = 1, n + 1 + thetaF(i) = thetac + (i - 1) * deltaf + end do CollarColatitudes + + ! Collars building + ypsi(1) = fpi*((sin(thetaf(2) / 2.0_dp))**2 - & + (sin(Thetaf(1) / 2.0_dp))**2) / vr + m(1) = int(ypsi(1) + 0.5_dp) + apsi(1) = ypsi(1) - m(1) + Collars: do i = 2, n + ypsi(i) = fpi * ((sin(thetaf(i+1) / 2.0_dp))**2 - & + (sin(thetaf(i ) / 2.0_dp))**2) / vr + m(i) = int(ypsi(i) + apsi(i-1) + 0.5_dp) + apsi(i) = ypsi(i) - m(i) + apsi(i-1) + end do Collars + + j = 1 + do i = 1, n + 1 + theta(i) = 2.0_dp * asin(sqrt((vr * j) / fpi)) + j = j + int(m(i) + 0.5_dp) + end do + + k = 2 + do i = 1, n + thetam = (theta(i) + theta(i+1)) / 2.0_dp + costm = cos(thetam) + sintm = sin(thetam) + TranslatePoints: do j = 1, int(m(i) + 0.5_dp) + omegaj = tpi * j / (1.0_dp * m(i)) + xtscor(nts+k) = x + r * sintm * cos(omegaj) + ytscor(nts+k) = y + r * sintm * sin(omegaj) + ztscor(nts+k) = z + r * costm + k=k+1 + end do TranslatePoints + end do + + end subroutine leopardi_grid + + subroutine write_visualization_file(nts, nesfp, xtscor, ytscor, ztscor, ar, xe, ye, ze, rin) + + !> writes a script for VMD + + !> Passed variables + integer(kind=regint_k), intent(in) :: nts, nesfp + real(kind=dp), intent(in) :: xtscor(nts), ytscor(nts), ztscor(nts), ar(nts), xe(nesfp), ye(nesfp), ze(nesfp), rin(nesfp) + !> Local variables + integer(kind=regint_k) :: file_unit + logical :: file_exists, file_opened + integer :: k + + file_unit = 221221_regint_k + file_exists = .false. + inquire(file = 'tsless_visual.vmd', exist = file_exists, opened = file_opened) + if (.not. file_opened) then + open(file_unit, & + file = 'tsless_visual.vmd', & + status = 'replace', & + form = 'formatted', & + access = 'sequential') + end if + + write(file_unit, *) 'mol load graphics name' + + do k = 1, nesfp + write(file_unit,1000) xe(k), ye(k), ze(k), rin(k) + enddo + + write(file_unit, *) 'graphics top color 3' + + do k = 1, nts + write(file_unit,1000) xtscor(k), ytscor(k), ztscor(k), 0.1_dp + end do + write(file_unit, *) + + close(file_unit) + +1000 format ('graphics top sphere {', 3F10.4, '} radius', F10.4, ' resolution 120') + + end subroutine write_visualization_file + + pure function cavity_surface(nts, ar) result(surface) + + !> Passed variables + integer(kind=regint_k), intent(in) :: nts + real(kind=dp), intent(in) :: ar(nts) + real(kind=dp) :: surface + !> Local variables + integer :: i + + surface = 0.0_dp + do i = 1, nts + surface = surface + ar(i) + end do + + end function cavity_surface + + pure function cavity_volume(nts, rsph, ar) result(volume) + + !> Passed variables + integer(kind=regint_k), intent(in) :: nts + real(kind=dp), intent(in) :: rsph(nts), ar(nts) + real(kind=dp) :: volume + !> Local variables + integer :: i + + volume = 0.0_dp + do i = 1, nts + volume = volume + ar(i) * rsph(i) + end do + volume = volume / 3.0_dp + + end function cavity_volume + +end module tsless_cavity diff --git a/src/tsless/tsless_lapack.f90 b/src/tsless/tsless_lapack.f90 new file mode 100644 index 000000000..aed4731bb --- /dev/null +++ b/src/tsless/tsless_lapack.f90 @@ -0,0 +1,158 @@ +! +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. +! +! This file is part of PCMSolver. +! +! PCMSolver is free software: you can redistribute it and/or modify +! it under the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! PCMSolver is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with PCMSolver. If not, see . +! +! For information on the complete list of contributors to the +! PCMSolver API, see: +! + +module tsless_lapack + +use tsless_precision + +implicit none + +public solve_linear_system + +private + +contains + + !> \brief Solve linear system of equations + !> \author Roberto Di Remigio + !> \year 2015 + !> \param[in] ndim dimension of the problem + !> \param[in] A system matrix + !> \param[in] b right-hand side + !> \param[out] x solution vector + !> Computes the explicit inverse by LU decomposition and then + !> performs matrix-vector multiplication A^-1 b = x + pure subroutine solve_linear_system(ndim, A, b, x) + + !> Passed variables + integer(kind=regint_k), intent(in) :: ndim + real(kind=dp), intent(inout) :: A(ndim, ndim) + real(kind=dp), intent(in) :: b(ndim) + real(kind=dp), intent(out) :: x(ndim) + !> Local variables + real(kind=dp), allocatable :: Ainverse(:, :) + integer(kind=regint_k) :: i, j + + !> Compute inverse by LU decomposition + allocate(Ainverse(ndim, ndim)); Ainverse = 0.0_dp + call inverse_LU(ndim, A, Ainverse) + + !> Solve linear system + x = 0.0_dp + do i = 1_regint_k, ndim + do j = 1_regint_k, ndim + x(i) = x(i) + Ainverse(i, j) * b(j) + end do + end do + + deallocate(Ainverse) + + end subroutine solve_linear_system + + !> \brief Computes inverse matrix + !> \author Alexander L. Godunov + !> \year 2009 + !> \param[in, out] A matrix to be inverted + !> \param[in, out] C inverse matrix + !> \param[in] ndim dimension of the matrix + !> The method is based on Doolittle LU factorization for Ax = b + !> notice that the contents of the original matrix will be destroyed + pure subroutine inverse_LU(ndim, A, C) + + !> Passed variables + integer(kind=regint_k), intent(in) :: ndim + real(kind=dp), intent(inout) :: A(ndim, ndim) + real(kind=dp), intent(inout) :: C(ndim, ndim) + !> Local variables + real(kind=dp), allocatable :: L(:, :), U(:, :) + real(kind=dp), allocatable :: b(:), d(:), x(:) + real(kind=dp) :: coeff + integer(kind=regint_k) :: i, j, k + + ! step 0: initialization for matrices L and U and b + allocate(L(ndim, ndim)); L = 0.0_dp + allocate(U(ndim, ndim)); U = 0.0_dp + allocate(b(ndim)); b = 0.0_dp + allocate(d(ndim)); d = 0.0_dp + allocate(x(ndim)); x = 0.0_dp + + ! step 1: forward elimination + do k = 1, ndim-1 + do i = k + 1, ndim + coeff = A(i, k) / A(k, k) + L(i, k) = coeff + do j = k + 1, ndim + A(i, j) = A(i, j) - coeff * A(k, j) + end do + end do + end do + + ! Step 2: prepare L and U matrices + ! L matrix is a matrix of the elimination coefficient + ! + the diagonal elements are 1.0 + do i = 1, ndim + L(i, i) = 1.0_dp + end do + ! U matrix is the upper triangular part of A + do j = 1, ndim + do i = 1, j + U(i, j) = A(i, j) + end do + end do + + ! Step 3: compute columns of the inverse matrix C + do k = 1, ndim + b(k) = 1.0_dp + d(1) = b(1) + ! Step 3a: Solve Ld=b using the forward substitution + do i = 2, ndim + d(i) = b(i) + do j = 1, i - 1 + d(i) = d(i) - L(i, j) * d(j) + end do + end do + ! Step 3b: Solve Ux=d using the back substitution + x(ndim) = d(ndim) / U(ndim, ndim) + do i = ndim - 1, 1, -1 + x(i) = d(i) + do j = ndim, i + 1, -1 + x(i) = x(i) - U(i, j) * x(j) + end do + x(i) = x(i) / U(i, i) + end do + ! Step 3c: fill the solutions x(ndim) into column k of C + do i = 1, ndim + C(i, k) = x(i) + end do + b(k) = 0.0_dp + end do + + deallocate(L) + deallocate(U) + deallocate(b) + deallocate(d) + deallocate(x) + + end subroutine inverse_LU + +end module tsless_lapack diff --git a/src/tsless/tsless_precision.f90 b/src/tsless/tsless_precision.f90 new file mode 100644 index 000000000..627605a45 --- /dev/null +++ b/src/tsless/tsless_precision.f90 @@ -0,0 +1,42 @@ +! +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. +! +! This file is part of PCMSolver. +! +! PCMSolver is free software: you can redistribute it and/or modify +! it under the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! PCMSolver is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with PCMSolver. If not, see . +! +! For information on the complete list of contributors to the +! PCMSolver API, see: +! + +module tsless_precision +! Read this: http://stackoverflow.com/a/3170438/2528668 +! and this: http://stackoverflow.com/a/3204981/2528668 + +implicit none + +! Integer types +! 32-bit integers +integer, parameter :: regint_k = selected_int_kind(8) +! 64-bit integers +integer, parameter :: largeint_k = selected_int_kind(18) + +! Real types +! Single-precision real +integer, parameter :: sp = kind(1.0) +! Double-precision real +integer, parameter :: dp = selected_real_kind(2*precision(1.0_sp)) + +end module tsless_precision diff --git a/src/tsless/tsless_symmetry.f90 b/src/tsless/tsless_symmetry.f90 new file mode 100644 index 000000000..92e1c2cb0 --- /dev/null +++ b/src/tsless/tsless_symmetry.f90 @@ -0,0 +1,354 @@ +! +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. +! +! This file is part of PCMSolver. +! +! PCMSolver is free software: you can redistribute it and/or modify +! it under the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! PCMSolver is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with PCMSolver. If not, see . +! +! For information on the complete list of contributors to the +! PCMSolver API, see: +! + +module tsless_symmetry +! NOTE: the ibtfun module is gone, as we can safely use Fortran +! standard intrinsic functions. +! Mapping between intrinsics and ibtfun: +! ibtand(i, j) <-> iand(i, j) +! ibtor(i, j) <-> ior(i, j) +! ibtshl(i, j) <-> ishft(i, j) +! ibtshr(i, j) <-> ishft(i, -j) !WARNING! +! ibtxor(i, j) <-> ieor(i, j) + +use, intrinsic :: iso_c_binding +use tsless_precision + +implicit none + +public get_pt +public build_point_group +! MINI-MANUAL +! A. Indexing of symmetry operations and their mapping to a bitstring: +! zyx Parity +! 0 000 E 1.0 +! 1 001 Oyz -1.0 +! 2 010 Oxz -1.0 +! 3 011 C2z 1.0 +! 4 100 Oxy -1.0 +! 5 101 C2y 1.0 +! 6 110 C2x 1.0 +! 7 111 i -1.0 +! B. Indexing of isymax matrix +! The isymax array contains the irrep to which the linear functions +! (first column) and the rotations (second column) belong. The indexing +! is given above. +! C. Indexing of the jsop array +! The jsop array contains the position at which the operation +! i (index given in point A) appears +! + +type, public :: point_group +! A type containing all you need to know about the point group + + ! String with the group name: + ! C1, C2, Cs, Ci, D2, C2v, C2h, D2h + character(len=3) :: group_name + ! Integer representing the group + ! 0, 1, 2, 3, 4, 5, 6, 7 + integer(kind=regint_k) :: group_int + ! Number of generators + integer(kind=regint_k) :: nr_generators + ! Number of not-trivial symmetry operations (2**nr_generators - 1) + integer(kind=regint_k) :: maxrep + ! group%isymax(i, 1): behaviour of principal axes under basic operations + ! (x-y-z) + ! group%isymax(i, 2): behaviour of principal rotations under basic operations + ! (Rx-Ry-Rz) + integer(kind=regint_k) :: isymax(3, 2) + ! Symmetry operations in the Abelian groups. + ! Bitstring: 1 coordinate changes sign under operation; + ! 0 coordinate does not change sign. + ! Of course, that's also the binary representation of + ! numbers from 0 to 7! + integer(kind=regint_k) :: jsop(0:7) + ! + integer(kind=regint_k) :: nr_rotations + ! + integer(kind=regint_k) :: nr_reflections + ! + integer(kind=regint_k) :: nr_inversion +end type point_group + +private +contains + + !> \brief returns parity of a bitstring + !> \param[in] bit_rep bitstring represenation of symmetry operator + !> PT is the parity of a bitstring: + !> 1 for an even number of ones: 000,011,110,101 + !> -1 for an odd number of ones: 001,010,100,111 + real(kind=dp) function get_pt(bit_rep) + + integer(kind=regint_k), intent(in) :: bit_rep + + real(kind=dp) :: pt(0:7) + + pt(0) = 1.0d0 + pt(1) = -1.0d0 + pt(2) = -1.0d0 + pt(3) = 1.0d0 + pt(4) = -1.0d0 + pt(5) = 1.0d0 + pt(6) = 1.0d0 + pt(7) = -1.0d0 + + get_pt = pt(bit_rep) + + end function get_pt + + !> \brief Builds point group given the generators + !> \param[in] nr_gen number of generators + !> \param[in] gen1 first generator + !> \param[in] gen2 second generator + !> \param[in] gen3 third generator + !> \param[in] print_unit logical unit for printing + !> Originally written by Trond Saue for DALTON/DIRAC + !> Copied and adapted by Roberto Di Remigio + !> The generators are encoded as bitstrings + function build_point_group(nr_gen, gen1, gen2, gen3, printer) result(pg) + + use, intrinsic :: iso_fortran_env, only: output_unit + + !> Passed variables + integer(kind=regint_k), intent(in) :: nr_gen + integer(kind=regint_k), intent(in) :: gen1, gen2, gen3 + integer, optional, intent(in) :: printer + !> Output variables + type(point_group) :: pg + !> Local variables + integer :: print_out + integer(kind=regint_k) :: maxrep + integer(kind=regint_k) :: isymax(3, 2) + integer(kind=regint_k) :: igen(3) + !> Integer representation of the rotations bitmaps + integer(kind=regint_k), parameter :: irots(3) = [3, 5, 6] + integer(kind=regint_k), parameter :: rots(3) = [6, 5, 3] + !> Integer representation of the reflections bitmaps + integer(kind=regint_k), parameter :: irefl(3) = [4, 2, 1] + !> Parity of the symmetry operations bitmaps + integer(kind=regint_k), parameter :: jpar(0:7) = [1, -1, -1, 1, -1, 1, 1, -1] + integer(kind=regint_k) :: i, j, k, l, i0, i1, i2, ind, ipos, bitmap + integer(kind=regint_k) :: nrots, nrefl, ninvc, igroup + integer(kind=regint_k) :: char_tab(0:7, 0:7) + logical :: lsymop(0:7) + integer(kind=regint_k) :: jsop(0:7), ipar(0:7) + character(3) :: group, rep(0:7) + character(3), parameter :: groups(0:7) = ['C1 ', 'C2 ', 'Cs ', & + 'Ci ', 'D2 ', 'C2v', & + 'C2h', 'D2h'] + character(3), parameter :: symop(0:7) = [' E ', 'Oyz', 'Oxz', & + 'C2z', 'Oxy', 'C2y', & + 'C2x', ' i '] + + if (present(printer)) then + print_out = printer + else + print_out = output_unit + end if + + isymax = 0 + igen = 0 + maxrep = 2**nr_gen - 1 + ! igen contains the bitmap for the generators + igen = [gen1, gen2, gen3] + ! Build isymax(:, 1) + ! determine to which irrep the translations belong to + ! Loop over Cartesian axes + do i = 1, 3 + bitmap = 0 + ! Loop over generators + do j = 1, nr_gen + ! Apply generators on Cartesian axes rots(i) and check the character + if (nint(get_pt(ior(igen(j), rots(i)))) == -1) then + ! Set the bitmap + bitmap = ibset(bitmap, j) + end if + end do + ! Right-shift the bitmap and assign to isymax + isymax(i, 1) = ishft(bitmap, -1) + end do + + ! Build isymax(:, 2) + ! determine to which irrep the rotations belong to + ! R_x = (y XOR z) and cyclic permutations + isymax(1, 2) = ieor(isymax(2, 1), isymax(3, 1)) + isymax(2, 2) = ieor(isymax(3, 1), isymax(1, 1)) + isymax(3, 2) = ieor(isymax(1, 1), isymax(2, 1)) + + ! Build the character table + lsymop = .false. + ! Activate all symmetry operations of the group + lsymop(0) = .true. + jsop(0) = 0 + ipar(0) = 1 + do i = 1, maxrep + i0 = iand(1_regint_k, i) * igen(1) + i1 = iand(1_regint_k, ishft(i, -1)) * igen(2) + i2 = iand(1_regint_k, ishft(i, -2)) * igen(3) + ind = ieor(ieor(i0, i1),i2) + lsymop(ind) = .true. + ipar(i) = jpar(ind) + end do + ! List group operations in preferred order + ! Identity, E + ind = 0 + jsop(ind) = 0 + ! Rotations + nrots = 0 + do i = 1, 3 + if (lsymop(irots(i))) then + ind = ind + 1 + jsop(ind) = irots(i) + nrots = nrots + 1 + end if + end do + ! Inversion + ninvc = 0 + if (lsymop(7)) then + ind = ind + 1 + jsop(ind) = 7 + ninvc = 1 + end if + ! Reflections + nrefl = 0 + do i = 1, 3 + if (lsymop(irefl(i))) then + ind = ind + 1 + jsop(ind) = irefl(i) + nrefl = nrefl + 1 + end if + end do + ! Classify group + ! ============== + ! tsaue - Here I have devised a highly empirical formula, but it works !!! + igroup = min(7, nint((4 * nrots + 8 * ninvc + 6 * nrefl) / 3.0)) + group = groups(igroup) + char_tab = 0 + ! Now generate the character table + do i = 0, maxrep + ! The character of the identity is always +1 + char_tab(0, i) = 1 + do j = 1, nr_gen + char_tab(igen(j), i) = nint(get_pt(iand(ishft(i,-(j-1)), 1_regint_k))) + do k = 1, (j-1) + ind = ieor(igen(j), igen(k)) + char_tab(ind, i) = char_tab(igen(j), i) * char_tab(igen(k), i) + do l = 1, (k-1) + char_tab(ieor(ind, igen(l)), i) = char_tab(ind, i) * char_tab(igen(l), i) + end do + end do + end do + end do + ! Classify irrep + do i = 0, maxrep + rep(i) = 'A ' + ipos = 2 + ! 1. Rotational symmetry + if (nrots == 3) then + ind = (1 - char_tab(jsop(1), i)) + (1 - char_tab(jsop(2), i))/2 + if (ind /= 0) then + rep(i)(1:1) = 'B' + rep(i)(2:2) = char(ichar('0') + ind) + ipos = 3 + end if + else if (nrots == 1) then + if (char_tab(jsop(1), i) == -1) then + rep(i)(1:1) = 'B' + end if + if (nrefl == 2) then + if (iand(ishft(jsop(1), -1), 1_regint_k) == 1) then + ind = 2 + else + ind = 3 + end if + if (char_tab(jsop(ind), i) == 1) then + rep(i)(2:2) = '1' + else + rep(i)(2:2) = '2' + end if + end if + else if (nrefl == 1) then + ! 2. Mirror symmetry + if (char_tab(jsop(1), i) == 1) then + rep(i)(2:2) = '''' + else if (char_tab(jsop(1), i) == -1) then + rep(i)(2:2) = '"' + end if + end if + ! 3. Inversion symmetry + if (ninvc == 1) then + ind = nrots + 1 + if (char_tab(jsop(ind), i) == 1) then + rep(i)(ipos:ipos) = 'g' + else + rep(i)(ipos:ipos) = 'u' + end if + end if + end do + ! Output + ! 1. Group name and generators + write(print_out, '(a, a3)') 'Point group: ', group + if (nr_gen > 0) then + write(print_out, '(/3x, a/)') '* The point group was generated by:' + do i = 1, nr_gen + if (symop(igen(i))(1:1) == 'C') then + write(print_out, '(6x, 3a)') 'Rotation about the ', symop(igen(i))(3:3),'-axis' + else if (symop(igen(i))(1:1) == 'O') then + write(print_out, '(6x, 3a)') 'Reflection in the ', symop(igen(i))(2:3),'-plane' + else + write(print_out, '(6x, a)') 'Inversion center' + end if + end do + ! 2. Group multiplication table + write(print_out,'(/3x, a/)') '* Group multiplication table' + write(print_out,'(8x, a1, 8(1x, a3, 1x))') '|', (symop(jsop(i)), i = 0, maxrep) + write(print_out,'(3x,a6,8a5)') '-----+', ('-----', i = 0, maxrep) + do i = 0, maxrep + write(print_out,'(4x, a3, 1x, a1, 8(1x, a3, 1x))') symop(jsop(i)), '|', & + (symop(ieor(jsop(i), jsop(j))), j = 0, maxrep) + end do + ! 3. Character table + write(print_out,'(/3x, a/)') '* Character table' + write(print_out,'(8x, a1, 8(1x, a3, 1x))') '|', (symop(jsop(j)), j = 0, maxrep) + write(print_out,'(3x, a6, 8a5)') '-----+', ('-----', i = 0, maxrep) + do i = 0, maxrep + write(print_out,'(4x, a3, 1x, a1, 8(1x, i3, 1x))') rep(i), '|', (char_tab(jsop(j), i), j=0, maxrep) + end do + ! 4. Direct product table + write(print_out,'(/3x, a/)') '* Direct product table' + write(print_out,'(8x, a1, 8(1x, a3, 1x))') '|', (rep(i), i = 0, maxrep) + write(print_out,'(3x, a6, 8a5)') '-----+', ('-----', i = 0, maxrep) + do i = 0, maxrep + write(print_out,'(3x, 1x, a3, 1x, a1, 8(1x, a3, 1x))') rep(i), '|', (rep(ieor(i, j)), j = 0, maxrep) + end do + end if + ! Fields: group name, group integer(kind=regint_k), number of generators, + ! number of nontrivial operations, isymax, jsop, + ! number of rotations, number of reflections, + ! number of inversions. + pg = point_group(group, igroup, nr_gen, maxrep, isymax, jsop, nrots, nrefl, ninvc) + + end function build_point_group + + end module tsless_symmetry diff --git a/src/tsless/tsless_weight_function.f90 b/src/tsless/tsless_weight_function.f90 new file mode 100644 index 000000000..a84663eb5 --- /dev/null +++ b/src/tsless/tsless_weight_function.f90 @@ -0,0 +1,248 @@ +! +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. +! +! This file is part of PCMSolver. +! +! PCMSolver is free software: you can redistribute it and/or modify +! it under the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! PCMSolver is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with PCMSolver. If not, see . +! +! For information on the complete list of contributors to the +! PCMSolver API, see: +! + +module tsless_weight_function + +use, intrinsic :: iso_c_binding +use tsless_precision + +implicit none + +!> \brief Describes a weight function +!> Data fiels have dimensions 2*nr_derivative+3 +type, public :: weight_function + !> Number of continuous derivatives + integer(kind=regint_k) :: nr_derivative + !> + real(kind=dp), pointer :: weight_1(:) + !> + real(kind=dp), pointer :: weight_2(:) + !> + real(kind=dp), pointer :: weight_3(:) +end type weight_function + +public create_weight_function +public destroy_weight_function +public scaling_factor + +private + +contains + + function create_weight_function(nderiv, normalization, dmin, printer) result(wfun) + + use, intrinsic :: iso_fortran_env, only: output_unit + + !> Passed variables + integer(kind=regint_k), intent(in) :: nderiv + integer(kind=regint_k), intent(in) :: normalization + real(kind=dp), intent(in) :: dmin + integer, optional, intent(in) :: printer + !> Output variables + type(weight_function) :: wfun + !> Local variables + integer :: print_out + + if (present(printer)) then + print_out = printer + else + print_out = output_unit + end if + + call allocate_weight_function(wfun, nderiv) + call compute(wfun, normalization, dmin, print_out) + + end function create_weight_function + + pure subroutine allocate_weight_function(wfun, nderiv) + + !> Passed variables + type(weight_function), intent(inout) :: wfun + integer(kind=regint_k), intent(in) :: nderiv + !> Local variables + integer(kind=regint_k) :: nparms + + wfun%nr_derivative = nderiv + nparms = 2_regint_k * nderiv + 3_regint_k + allocate(wfun%weight_1(nparms)); wfun%weight_1 = 0.0_dp + allocate(wfun%weight_2(nparms)); wfun%weight_2 = 0.0_dp + allocate(wfun%weight_3(nparms)); wfun%weight_3 = 0.0_dp + + end subroutine allocate_weight_function + + pure subroutine destroy_weight_function(wfun) + + !> Passed variables + type(weight_function), intent(inout) :: wfun + + wfun%nr_derivative = 0_regint_k + deallocate(wfun%weight_1) + deallocate(wfun%weight_2) + deallocate(wfun%weight_3) + + end subroutine destroy_weight_function + + !> \brief + !> \author Christian S. Pomelli + !> \param[inout] wfun weight function + !> \param[in] nord maximum number of continuous derivatives + !> \param[in] ifun + !> \param[in] dmin minimal distance between sampling points + !> \param[in] printer logical unit for printing + !> ifun = 0 non-normalized function (gamma) + !> ifun = 1 normalized function (omega) + subroutine compute(wfun, ifun, dmin, printer) + + use, intrinsic :: iso_fortran_env, only: output_unit + + use tsless_lapack, only: solve_linear_system + + !> Passed variables + type(weight_function), intent(inout) :: wfun + integer(kind=regint_k), intent(in) :: ifun + real(kind=dp), intent(in) :: dmin + integer, optional, intent(in) :: printer + !> Local variables + real(kind=dp), allocatable :: C(:, :) + real(kind=dp), allocatable :: b(:) + integer :: i, j, k, n, m + integer :: print_out + + if (present(printer)) then + print_out = printer + else + print_out = output_unit + end if + + ! check parameters + if ((ifun .lt. 0) .or. (ifun .gt. 1)) Then + write(print_out, *) 'Unknown IFUN' + stop + end if + + ! m = number of parameters + n = wfun%nr_derivative + if (ifun .eq. 0) then + m = 2*n + 2 + else + m = 2*n + 3 + end if + !> Allocate, clean up and fill coefficient matrix C + allocate(C(m, m)); C = 0.0_dp + ! l = 1 boundary conditions + do j = 0_regint_k, n + do k = j, m - 1_regint_k + C(j+1, k+1) = coefficients(j, k) + end do + end do + ! l = d boundary conditions + do j = 0_regint_k, n + do k = j, m - 1_regint_k + C(j+2+n, k+1) = coefficients(j, k) + C(j+2+n, k+1) = C(j+2+n, k+1) * dmin**(k-j) + end do + end do + ! normalization condition + if (ifun .eq. 1_regint_k) then + do k = 0_regint_k, m - 1_regint_k + C(m, k+1) = (1.0_dp - dmin**(k+1)) / (1.0_dp * (k+1)) + end do + end if + ! Allocate, clean up and fill b vector (RHS) + allocate(b(m)); b = 0.0_dp + b(1) = 1.0_dp + if (ifun .eq. 1) b(m) = 1.0_dp - dmin + ! 2. Solve linear system of equations + call solve_linear_system(ndim=m, A=C, b=b, x=wfun%weight_1) + ! write function and derivatives + do i = 1_regint_k, m - 1_regint_k + wfun%weight_2(i) = i * wfun%weight_1(i+1) + end do + + do i = 1_regint_k, m - 2_regint_k + wfun%weight_3(i) = i * wfun%weight_2(i+1) + end do + + write(print_out, *) 'ifun = ', ifun + do i = 1_regint_k, m + write(print_out, '(a, i2, a, 3F15.4)') 'c(', i, ') = ', wfun%weight_1(i), wfun%weight_2(i), wfun%weight_3(i) + end do + write(print_out, *) '-------------------' + + deallocate(b) + deallocate(C) + + end subroutine compute + + !> \brief Calculates scaling factor for weight of point + !> \author Christian S. Pomelli + !> \param[in] wfun weight function + !> \param[in] arg argument of the regularized step function + !> \param[in] dmin minimal distance between sampling points + pure function scaling_factor(wfun, arg, dmin) result(scaling) + + !> Input variables + type(weight_function), intent(in) :: wfun + real(kind=dp), intent(in) :: arg + real(kind=dp), intent(in) :: dmin + !> Output variables + real(kind=dp) :: scaling(3) + !> Local variables + integer :: i, nr_params + real(kind=dp) :: tn + + ! Simply a polynomial. later add derivatives + scaling = 0.0_dp + tn = 1.0_dp + ! number of polynomial parameters + nr_params = 2_regint_k * wfun%nr_derivative + 3_regint_k + do i = 1_regint_k, nr_params + scaling(1) = scaling(1) + wfun%weight_1(i) * tn + scaling(2) = scaling(2) + wfun%weight_2(i) * tn + scaling(3) = scaling(3) + wfun%weight_3(i) * tn + tn = tn * arg + end do + + end function scaling_factor + + pure function coefficients(j, k) result(f) + + !> Passed variables + integer(kind=regint_k), intent(in) :: j, k + real(kind=dp) :: f + !> Local variables + integer :: i + real(kind=dp) :: dkk + + f = 1.0_dp + dkk = dble(k) + if (j .ne. 0) then + do i = 1, j + f = f * dkk + dkk = dkk - 1.0_dp + end do + end if + + end function coefficients + +end module tsless_weight_function diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3982e576c..897932968 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -11,16 +11,17 @@ include_directories( ${PROJECT_SOURCE_DIR}/tests ) -add_subdirectory(input) -add_subdirectory(numerical_quadrature) -add_subdirectory(gepol) -add_subdirectory(dielectric_profile) -add_subdirectory(green) +add_subdirectory(benchmark) add_subdirectory(bi_operators) add_subdirectory(cpcm) +add_subdirectory(dielectric_profile) +add_subdirectory(gepol) +add_subdirectory(green) add_subdirectory(iefpcm) +add_subdirectory(input) +add_subdirectory(numerical_quadrature) +add_subdirectory(tsless) add_subdirectory(utils) -add_subdirectory(benchmark) add_executable(unit_tests unit_tests.cpp $ @@ -31,6 +32,7 @@ add_executable(unit_tests unit_tests.cpp $ $ $ + $ $ ) if(BUILD_CUSTOM_BOOST) diff --git a/tests/TestingMolecules.hpp b/tests/TestingMolecules.hpp index ce51f7d3e..5ed18726e 100644 --- a/tests/TestingMolecules.hpp +++ b/tests/TestingMolecules.hpp @@ -35,6 +35,10 @@ #include "utils/Symmetry.hpp" namespace pcm { +/*! Returns the hydrogen fluoride molecule + */ +inline Molecule HF(); + /*! Returns the ammonia molecule */ inline Molecule NH3(); @@ -123,6 +127,34 @@ template Molecule dummy(double radius, const Eigen::Vector3d & cente return dummy; }; +Molecule HF() { + int nAtoms = 2; + + Eigen::Vector3d F(0.000000000, 0.00000000, 0.08729478); + Eigen::Vector3d H(0.000000000, 0.00000000, -1.64558444); + + Eigen::MatrixXd geom(3, nAtoms); + geom.col(0) = F.transpose(); + geom.col(1) = H.transpose(); + Eigen::Vector2d charges, masses; + charges << 9.0, 1.0; + masses << 18.9984030, 1.0078250; + std::vector atoms; + atoms.push_back(Atom("Fluorine", "F", charges(0), masses(0), 2.777897403, F, 1.0)); + atoms.push_back(Atom("Hydrogen", "H", charges(1), masses(1), 2.267671349, H, 1.0)); + + std::vector spheres; + Sphere sph1(F, 2.777897403); + Sphere sph2(H, 2.267671349); + spheres.push_back(sph1); + spheres.push_back(sph2); + + // C1 + Symmetry pGroup = buildGroup(0, 0, 0, 0); + + return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup); +}; + Molecule NH3() { int nAtoms = 4; diff --git a/tests/cpcm/CMakeLists.txt b/tests/cpcm/CMakeLists.txt index c5bdbbdb3..90a26eb95 100644 --- a/tests/cpcm/CMakeLists.txt +++ b/tests/cpcm/CMakeLists.txt @@ -6,6 +6,8 @@ add_library(cpcm-tests OBJECT cpcm_gepol-point.cpp cpcm_gepol-point_from-file.cpp cpcm_symmetry.cpp + cpcm_tsless-point.cpp + cpcm_tsless-NH3.cpp ) if(BUILD_CUSTOM_BOOST) add_dependencies(cpcm-tests custom_boost) @@ -96,3 +98,22 @@ add_Catch_test( REFERENCE_FILES ASC-cpcm_gepol-C2H4_D2h.npy ) + +# FIXME Reference files for these two tests +# cpcm_tsless-point.cpp test +add_Catch_test( + NAME + cpcm_tsless-point + LABELS + cpcm + cpcm_tsless-point + ) + +# cpcm_tsless-NH3.cpp test +add_Catch_test( + NAME + cpcm_tsless-NH3 + LABELS + cpcm + cpcm_tsless-NH3 + ) diff --git a/tests/cpcm/cpcm_tsless-NH3.cpp b/tests/cpcm/cpcm_tsless-NH3.cpp new file mode 100644 index 000000000..56a5f9d5c --- /dev/null +++ b/tests/cpcm/cpcm_tsless-NH3.cpp @@ -0,0 +1,86 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PCMSolver is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include "Config.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "bi_operators/Collocation.hpp" +#include "cavity/TsLessCavity.hpp" +#include "green/DerivativeTypes.hpp" +#include "green/UniformDielectric.hpp" +#include "green/Vacuum.hpp" +#include "solver/CPCMSolver.hpp" +#include "utils/Molecule.hpp" +#include "utils/Symmetry.hpp" + +using namespace pcm; +using bi_operators::Collocation; +using cavity::TsLessCavity; +using green::UniformDielectric; +using green::Vacuum; +using solver::CPCMSolver; + +/*! \class CPCMSolver + * \test \b NH3TsLess tests CPCMSolver using ammonia and a TsLess cavity + */ +TEST_CASE("Test solver for the C-PCM with NH3 molecule and a TsLess cavity", + "[solver][cpcm][cpcm_tsless-NH3]") { + Molecule molec = NH3(); + + double area = 0.08; + double minDistance = 0.1; + double probeRadius = 0.0; + int derOrder = 8; + double minRadius = 100.0; + TsLessCavity cavity = + TsLessCavity(molec, area, probeRadius, minRadius, minDistance, derOrder); + + double permittivity = 78.39; + Vacuum<> gfInside; + UniformDielectric<> gfOutside(permittivity); + bool symm = true; + double correction = 0.0; + + Collocation S; + + CPCMSolver solver(symm, correction); + solver.buildSystemMatrix(cavity, gfInside, gfOutside, S); + + double Ncharge = 7.0; + double Hcharge = 1.0; + size_t size = cavity.size(); + Eigen::VectorXd fake_mep = computeMEP(molec, cavity.elements()); + // The total ASC for a conductor is -Q + // for CPCM it will be -Q*[(epsilon-1)/epsilon + correction] + Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size); + fake_asc = solver.computeCharge(fake_mep); + double totalASC = + -(Ncharge + 3.0 * Hcharge) * (permittivity - 1) / (permittivity + correction); + double totalFakeASC = fake_asc.sum(); + CAPTURE(totalASC - totalFakeASC); + REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03)); +} diff --git a/tests/cpcm/cpcm_tsless-point.cpp b/tests/cpcm/cpcm_tsless-point.cpp new file mode 100644 index 000000000..1a267e675 --- /dev/null +++ b/tests/cpcm/cpcm_tsless-point.cpp @@ -0,0 +1,140 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PCMSolver is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include "Config.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "bi_operators/Collocation.hpp" +#include "cavity/TsLessCavity.hpp" +#include "green/DerivativeTypes.hpp" +#include "green/UniformDielectric.hpp" +#include "green/Vacuum.hpp" +#include "solver/CPCMSolver.hpp" + +using namespace pcm; +using bi_operators::Collocation; +using cavity::TsLessCavity; +using green::UniformDielectric; +using green::Vacuum; +using solver::CPCMSolver; + +SCENARIO("Test solver for the C-PCM for a point charge and a TsLess cavity", + "[solver][cpcm][cpcm_tsless-point]") { + GIVEN("An isotropic environment modelled as a perfect conductor and a point " + "charge") { + double permittivity = 78.39; + Vacuum<> gfInside; + UniformDielectric<> gfOutside(permittivity); + bool symm = true; + double correction = 0.0; + + Collocation S; + + double charge = 8.0; + double totalASC = -charge * (permittivity - 1) / (permittivity + correction); + + /*! \class CPCMSolver + * \test \b pointChargeTsLess tests CPCMSolver using a point charge with a + * TsLess cavity The point charge is at the origin. + */ + WHEN("the point charge is located at the origin") { + Molecule point = dummy<0>(2.929075493); + double area = 0.4; + double probeRadius = 0.0; + double minRadius = 100.0; + double minDistance = 0.1; + int derOrder = 4; + TsLessCavity cavity = + TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder); + cavity.saveCavity("tsless_point.npz"); + + CPCMSolver solver(symm, correction); + solver.buildSystemMatrix(cavity, gfInside, gfOutside, S); + + size_t size = cavity.size(); + Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge); + for (size_t i = 0; i < size; ++i) { + INFO("fake_mep(" << i << ") = " << fake_mep(i)); + } + + THEN("the apparent surface charge is") { + Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size); + fake_asc = solver.computeCharge(fake_mep); + double totalFakeASC = fake_asc.sum(); + + for (size_t i = 0; i < size; ++i) { + INFO("fake_asc(" << i << ") = " << fake_asc(i)); + } + + CAPTURE(totalASC); + CAPTURE(totalFakeASC); + CAPTURE(totalASC - totalFakeASC); + REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03)); + } + } + + /*! \class CPCMSolver + * \test \b pointChargeShiftedTsLess tests CPCMSolver using a point charge with + * a TsLess cavity The point charge is away from the origin. + */ + AND_WHEN("the point charge is located away from the origin") { + Eigen::Vector3d origin = 100 * Eigen::Vector3d::Random(); + Molecule point = dummy<0>(2.929075493, origin); + double area = 0.4; + double probeRadius = 0.0; + double minRadius = 100.0; + double minDistance = 0.1; + int derOrder = 4; + TsLessCavity cavity = + TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder); + + CPCMSolver solver(symm, correction); + solver.buildSystemMatrix(cavity, gfInside, gfOutside, S); + + size_t size = cavity.size(); + Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge, origin); + for (size_t i = 0; i < size; ++i) { + INFO("fake_mep(" << i << ") = " << fake_mep(i)); + } + + THEN("the surface charge is") { + Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size); + fake_asc = solver.computeCharge(fake_mep); + double totalFakeASC = fake_asc.sum(); + + for (size_t i = 0; i < size; ++i) { + INFO("fake_asc(" << i << ") = " << fake_asc(i)); + } + + CAPTURE(totalASC); + CAPTURE(totalFakeASC); + CAPTURE(totalASC - totalFakeASC); + REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03)); + } + } + } +} diff --git a/tests/gepol/CMakeLists.txt b/tests/gepol/CMakeLists.txt index 7b139a155..0b1296eb7 100644 --- a/tests/gepol/CMakeLists.txt +++ b/tests/gepol/CMakeLists.txt @@ -11,6 +11,7 @@ add_library(gepol-tests OBJECT gepol_point.cpp gepol_point_from-file.cpp gepol_point_symmetry.cpp + gepol_HF.cpp ) if(BUILD_CUSTOM_BOOST) add_dependencies(gepol-tests custom_boost) @@ -129,3 +130,12 @@ add_Catch_test( DEPENDS gepol_NH3 ) + +# gepol_HF.cpp test +add_Catch_test( + NAME + gepol_HF + LABELS + gepol + gepol_HF + ) diff --git a/tests/gepol/gepol_HF.cpp b/tests/gepol/gepol_HF.cpp new file mode 100644 index 000000000..5919f1207 --- /dev/null +++ b/tests/gepol/gepol_HF.cpp @@ -0,0 +1,76 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PCMSolver is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "cavity/GePolCavity.hpp" +#include "utils/Molecule.hpp" +#include "utils/Symmetry.hpp" + +using namespace pcm; +using cavity::GePolCavity; + +TEST_CASE("GePol cavity for an hydrogen fluoride molecule", "[gepol][gepol_HF]") { + Molecule molec = HF(); + double area = 0.02 / bohr2ToAngstrom2(); + double probeRadius = 0.0 / bohrToAngstrom(); + double minRadius = 0.2 / bohrToAngstrom(); + GePolCavity cavity(molec, area, probeRadius, minRadius); + + /*! \class GePolCavity + * \test \b GePolCavityHFTest_size tests GePol cavity size for ammonia + */ + SECTION("Test size") { + size_t ref_size = 1688; + size_t size = cavity.size(); + REQUIRE(size == ref_size); + } + + /*! \class GePolCavity + * \test \b GePolCavityHFTest_area tests GePol cavity surface area for ammonia + */ + SECTION("Test surface area") { + double ref_area = 110.64517236179323; + double area = cavity.elementArea().sum(); + REQUIRE(area == Approx(ref_area)); + } + + /*! \class GePolCavity + * \test \b GePolCavityHFTest_volume tests GePol cavity volume for ammonia + */ + SECTION("Test volume") { + double ref_volume = 105.98002389543836; + Eigen::Matrix3Xd elementCenter = cavity.elementCenter(); + Eigen::Matrix3Xd elementNormal = cavity.elementNormal(); + double volume = 0; + for (size_t i = 0; i < cavity.size(); ++i) { + volume += + cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i)); + } + volume /= 3; + REQUIRE(volume == Approx(ref_volume)); + } +} diff --git a/tests/iefpcm/CMakeLists.txt b/tests/iefpcm/CMakeLists.txt index c672cb086..ec878837c 100644 --- a/tests/iefpcm/CMakeLists.txt +++ b/tests/iefpcm/CMakeLists.txt @@ -11,6 +11,8 @@ add_library(iefpcm-tests OBJECT iefpcm_gepol-point.cpp iefpcm_gepol-point_from-file.cpp iefpcm_symmetry.cpp + iefpcm_tsless-point.cpp + iefpcm_tsless-NH3.cpp ) if(BUILD_CUSTOM_BOOST) add_dependencies(iefpcm-tests custom_boost) @@ -151,3 +153,23 @@ add_Catch_test( DEPENDS iefpcm_gepol-point ) + +# iefpcm_tsless-point.cpp test +add_Catch_test( + NAME + iefpcm_tsless-point + LABELS + solver + iefpcm + iefpcm_tsless-point + ) + +# iefpcm_tsless-NH3.cpp test +add_Catch_test( + NAME + iefpcm_tsless-NH3 + LABELS + solver + iefpcm + iefpcm_tsless-NH3 + ) diff --git a/tests/iefpcm/iefpcm_tsless-NH3.cpp b/tests/iefpcm/iefpcm_tsless-NH3.cpp new file mode 100644 index 000000000..1340d2891 --- /dev/null +++ b/tests/iefpcm/iefpcm_tsless-NH3.cpp @@ -0,0 +1,83 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PCMSolver is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include "Config.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "bi_operators/Collocation.hpp" +#include "cavity/TsLessCavity.hpp" +#include "green/DerivativeTypes.hpp" +#include "green/UniformDielectric.hpp" +#include "green/Vacuum.hpp" +#include "solver/IEFSolver.hpp" +#include "utils/Molecule.hpp" +#include "utils/Symmetry.hpp" + +using namespace pcm; +using bi_operators::Collocation; +using cavity::TsLessCavity; +using green::UniformDielectric; +using green::Vacuum; +using solver::IEFSolver; + +/*! \class IEFSolver + * \test \b NH3TsLess tests IEFSolver using ammonia and a TsLess cavity + */ +TEST_CASE("Test solver for the IEFPCM with NH3 molecule and a TsLess cavity", + "[solver][iefpcm][iefpcm_tsless-NH3]") { + Molecule molec = NH3(); + + double area = 0.08; + double minDistance = 0.1; + double probeRadius = 0.0; + int derOrder = 8; + double minRadius = 100.0; + TsLessCavity cavity = + TsLessCavity(molec, area, probeRadius, minRadius, minDistance, derOrder); + + double permittivity = 78.39; + Vacuum<> gfInside; + UniformDielectric<> gfOutside(permittivity); + + Collocation op; + + bool symm = true; + IEFSolver solver(symm); + solver.buildSystemMatrix(cavity, gfInside, gfOutside, op); + + double Ncharge = 7.0; + double Hcharge = 1.0; + size_t size = cavity.size(); + Eigen::VectorXd fake_mep = computeMEP(molec, cavity.elements()); + // The total ASC for a dielectric is -Q*[(epsilon-1)/epsilon] + Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size); + fake_asc = solver.computeCharge(fake_mep); + double totalASC = -(Ncharge + 3.0 * Hcharge) * (permittivity - 1) / permittivity; + double totalFakeASC = fake_asc.sum(); + CAPTURE(totalASC - totalFakeASC); + REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03)); +} diff --git a/tests/iefpcm/iefpcm_tsless-point.cpp b/tests/iefpcm/iefpcm_tsless-point.cpp new file mode 100644 index 000000000..cc22424ad --- /dev/null +++ b/tests/iefpcm/iefpcm_tsless-point.cpp @@ -0,0 +1,138 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PCMSolver is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include "Config.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "bi_operators/Collocation.hpp" +#include "cavity/TsLessCavity.hpp" +#include "green/DerivativeTypes.hpp" +#include "green/UniformDielectric.hpp" +#include "green/Vacuum.hpp" +#include "solver/IEFSolver.hpp" + +using namespace pcm; +using bi_operators::Collocation; +using cavity::TsLessCavity; +using green::UniformDielectric; +using green::Vacuum; +using solver::IEFSolver; + +SCENARIO("Test solver for the IEFPCM for a point charge and a TsLess cavity", + "[solver][iefpcm][iefpcm_tsless-point]") { + GIVEN("An isotropic environment and a point charge") { + double permittivity = 78.39; + Vacuum<> gfInside; + UniformDielectric<> gfOutside(permittivity); + bool symm = true; + + Collocation op; + + double charge = 8.0; + double totalASC = -charge * (permittivity - 1) / permittivity; + + /*! \class IEFSolver + * \test \b pointChargeTsLess tests IEFSolver using a point charge with a TsLess + * cavity The point charge is at the origin. + */ + WHEN("the point charge is located at the origin") { + Molecule point = dummy<0>(2.929075493); + double area = 0.4; + double probeRadius = 0.0; + double minRadius = 100.0; + double minDistance = 0.1; + int derOrder = 4; + TsLessCavity cavity = + TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder); + cavity.saveCavity("tsless_point.npz"); + + IEFSolver solver(symm); + solver.buildSystemMatrix(cavity, gfInside, gfOutside, op); + + size_t size = cavity.size(); + Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge); + // The total ASC for a dielectric is -Q*[(epsilon-1)/epsilon] + Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size); + fake_asc = solver.computeCharge(fake_mep); + + for (size_t i = 0; i < size; ++i) { + INFO("fake_mep(" << i << ") = " << fake_mep(i)); + } + for (size_t i = 0; i < size; ++i) { + INFO("fake_asc(" << i << ") = " << fake_asc(i)); + } + + double totalFakeASC = fake_asc.sum(); + THEN("the apparent surface charge is") { + CAPTURE(totalASC); + CAPTURE(totalFakeASC); + CAPTURE(totalASC - totalFakeASC); + REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03)); + } + } + + /*! \class IEFSolver + * \test \b pointChargeShiftedTsLess tests IEFSolver using a point charge with a + * TsLess cavity The point charge is away from the origin. + */ + AND_WHEN("the point charge is located away from the origin") { + Eigen::Vector3d origin = 100 * Eigen::Vector3d::Random(); + Molecule point = dummy<0>(2.929075493, origin); + double area = 0.4; + double probeRadius = 0.0; + double minRadius = 100.0; + double minDistance = 0.1; + int derOrder = 4; + TsLessCavity cavity = + TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder); + + IEFSolver solver(symm); + solver.buildSystemMatrix(cavity, gfInside, gfOutside, op); + + size_t size = cavity.size(); + Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge, origin); + // The total ASC for a dielectric is -Q*[(epsilon-1)/epsilon] + Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size); + fake_asc = solver.computeCharge(fake_mep); + + for (size_t i = 0; i < size; ++i) { + INFO("fake_mep(" << i << ") = " << fake_mep(i)); + } + for (size_t i = 0; i < size; ++i) { + INFO("fake_asc(" << i << ") = " << fake_asc(i)); + } + + double totalFakeASC = fake_asc.sum(); + THEN("the surface charge is") { + CAPTURE(totalASC); + CAPTURE(totalFakeASC); + CAPTURE(totalASC - totalFakeASC); + REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03)); + } + } + } +} diff --git a/tests/tsless/CMakeLists.txt b/tests/tsless/CMakeLists.txt new file mode 100644 index 000000000..e536b8ae8 --- /dev/null +++ b/tests/tsless/CMakeLists.txt @@ -0,0 +1,36 @@ +add_library(tsless-tests OBJECT + tsless_HF.cpp + tsless_NH3.cpp + tsless_point.cpp + ) +if(BUILD_CUSTOM_BOOST) + add_dependencies(tsless-tests custom_boost) +endif() + +# tsless_HF.cpp test +add_Catch_test( + NAME + tsless_HF + LABELS + tsless + tsless_HF + ) + +# tsless_point.cpp test +add_Catch_test( + NAME + tsless_point + LABELS + tsless + tsless_point + ) + +# tsless_NH3.cpp test +add_Catch_test( + NAME + tsless_NH3 + LABELS + tsless + tsless_NH3 + ) + diff --git a/tests/tsless/tsless_HF.cpp b/tests/tsless/tsless_HF.cpp new file mode 100644 index 000000000..c0245f3ec --- /dev/null +++ b/tests/tsless/tsless_HF.cpp @@ -0,0 +1,76 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PCMSolver is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "cavity/TsLessCavity.hpp" + +using namespace pcm; +using cavity::TsLessCavity; + +TEST_CASE("TsLess cavity for an hydrogen fluoride molecule", "[tsless][tsless_HF]") { + Molecule molec = HF(); + double area = 0.02 / bohr2ToAngstrom2(); + double probeRadius = 0.0 / bohrToAngstrom(); + double minRadius = 0.2 / bohrToAngstrom(); + double minDistance = 0.01; + int derOrder = 4; + TsLessCavity cavity(molec, area, probeRadius, minRadius, minDistance, derOrder); + + /*! \class TsLessCavity + * \test \b TsLessCavityHFTest_size tests TsLess cavity size for ammonia + */ + SECTION("Test size") { + size_t ref_size = 1498; + size_t size = cavity.size(); + REQUIRE(size == ref_size); + } + + /*! \class TsLessCavity + * \test \b TsLessCavityHFTest_area tests TsLess cavity surface area for ammonia + */ + SECTION("Test surface area") { + double ref_area = 111.07794350824591; + double area = cavity.elementArea().sum(); + REQUIRE(area == Approx(ref_area)); + } + + /*! \class TsLessCavity + * \test \b TsLessCavityHFTest_volume tests TsLess cavity volume for ammonia + */ + SECTION("Test volume") { + double ref_volume = 106.59560252994531; + Eigen::Matrix3Xd elementCenter = cavity.elementCenter(); + Eigen::Matrix3Xd elementNormal = cavity.elementNormal(); + double volume = 0; + for (size_t i = 0; i < cavity.size(); ++i) { + volume += + cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i)); + } + volume /= 3; + REQUIRE(volume == Approx(ref_volume)); + } +} diff --git a/tests/tsless/tsless_NH3.cpp b/tests/tsless/tsless_NH3.cpp new file mode 100644 index 000000000..ec611bb2c --- /dev/null +++ b/tests/tsless/tsless_NH3.cpp @@ -0,0 +1,77 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PCMSolver is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "cavity/TsLessCavity.hpp" + +using namespace pcm; +using cavity::TsLessCavity; + +TEST_CASE("TsLess cavity for an ammonia molecule", "[tsless][tsless_NH3]") { + Molecule molec = NH3(); + double area = 0.02 / bohr2ToAngstrom2(); + double probeRadius = 0.0 / bohrToAngstrom(); + double minRadius = 0.2 / bohrToAngstrom(); + double minDistance = 0.01; + int derOrder = 4; + TsLessCavity cavity(molec, area, probeRadius, minRadius, minDistance, derOrder); + cavity.saveCavity("nh3.npz"); + + /*! \class TsLessCavity + * \test \b TsLessCavityNH3Test_size tests TsLess cavity size for ammonia + */ + SECTION("Test size") { + size_t ref_size = 2042; + size_t size = cavity.size(); + REQUIRE(size == ref_size); + } + + /*! \class TsLessCavity + * \test \b TsLessCavityNH3Test_area tests TsLess cavity surface area for ammonia + */ + SECTION("Test surface area") { + double ref_area = 147.59900736316496; + double area = cavity.elementArea().sum(); + REQUIRE(area == Approx(ref_area)); + } + + /*! \class TsLessCavity + * \test \b TsLessCavityNH3Test_volume tests TsLess cavity volume for ammonia + */ + SECTION("Test volume") { + double ref_volume = 153.7079082474709; + Eigen::Matrix3Xd elementCenter = cavity.elementCenter(); + Eigen::Matrix3Xd elementNormal = cavity.elementNormal(); + double volume = 0; + for (size_t i = 0; i < cavity.size(); ++i) { + volume += + cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i)); + } + volume /= 3; + REQUIRE(volume == Approx(ref_volume)); + } +} diff --git a/tests/tsless/tsless_point.cpp b/tests/tsless/tsless_point.cpp new file mode 100644 index 000000000..6d032e171 --- /dev/null +++ b/tests/tsless/tsless_point.cpp @@ -0,0 +1,76 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PCMSolver is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "cavity/TsLessCavity.hpp" + +using namespace pcm; +using cavity::TsLessCavity; + +TEST_CASE("TsLess cavity for a single sphere", "[tsless][tsless_point]") { + Molecule point = dummy<0>(); + double area = 0.4; + double minDistance = 0.1; + double probeRadius = 0.0; + int derOrder = 4; + TsLessCavity cavity(point, area, probeRadius, 100.0, minDistance, derOrder); + + /*! \class TsLessCavity + * \test \b TsLessCavityTest_size tests TsLess cavity size for a point charge + */ + SECTION("Test size") { + size_t ref_size = 32; + size_t size = cavity.size(); + REQUIRE(ref_size == size); + } + + /*! \class TsLessCavity + * \test \b TsLessCavityTest_area tests TsLess cavity surface area for a point + * charge + */ + SECTION("Test surface area") { + double ref_area = 4.0 * M_PI * pow(1.0, 2); + double area = cavity.elementArea().sum(); + REQUIRE(ref_area == Approx(area)); + } + + /*! \class TsLessCavity + * \test \b TsLessCavityTest_volume tests TsLess cavity volume for a point charge + */ + SECTION("Test volume") { + double ref_volume = 4.0 * M_PI * pow(1.0, 3) / 3.0; + Eigen::Matrix3Xd elementCenter = cavity.elementCenter(); + Eigen::Matrix3Xd elementNormal = cavity.elementNormal(); + double volume = 0; + for (size_t i = 0; i < cavity.size(); ++i) { + volume += + cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i)); + } + volume /= 3; + REQUIRE(ref_volume == Approx(volume)); + } +} diff --git a/tools/pcmparser.py b/tools/pcmparser.py index 268020cea..ef7f4c2ab 100644 --- a/tools/pcmparser.py +++ b/tools/pcmparser.py @@ -168,47 +168,57 @@ def setup_keywords(): # Cavity section cavity = Section('CAVITY', callback=verify_cavity) # Type of the cavity - # Valid values: GEPOL and RESTART + # Valid values: GEPOL, TSLESS and RESTART cavity.add_kw('TYPE', 'STR') # Name of the file containing data for restarting a cavity # Valid for: Restart cavity # Default: empty string cavity.add_kw('NPZFILE', 'STR', '') # Average area (in au^2) - # Valid for: GePol + # Valid for: GePol and TsLess cavities # Valid values: double strictly greater than 0.01 au^2 # Default: 0.3 au^2 cavity.add_kw('AREA', 'DBL', 0.3) + # Minimal distance between sampling points (in au) + # Valid for: TsLess cavity + # Valid values: double + # Default: 0.1 au + cavity.add_kw('MINDISTANCE', 'DBL', 0.1) + # Derivative order for the weight function + # Valid for: TsLess cavity + # Valid values: integer + # Default: 4 + cavity.add_kw('DERORDER', 'INT', 4) # Scaling of the atomic radii - # Valid for: GePol + # Valid for: GePol and TsLess # Valid values: boolean # Default: True cavity.add_kw('SCALING', 'BOOL', True) # Built-in radii set - # Valid for: GePol + # Valid for: GePol and TsLess # Valid values: BONDI, UFF or ALLINGER # Default: BONDI cavity.add_kw('RADIISET', 'STR', 'BONDI') # Minimal radius of added spheres (in au) - # Valid for: GePol + # Valid for: GePol and TsLess # Valid values: double greater than 0.4 au # Default: 100.0 au (no added spheres) cavity.add_kw('MINRADIUS', 'DBL', 100.0) # Spheres geometry creation mode - # Valid for: GePol + # Valid for: GePol and TsLess # Valid values: EXPLICIT, ATOMS or IMPLICIT # Default: IMPLICIT cavity.add_kw('MODE', 'STR', 'IMPLICIT') # List of atoms with custom radius - # Valid for: GePol in MODE=ATOMS + # Valid for: GePol and TsLess in MODE=ATOMS # Valid values: array of integers cavity.add_kw('ATOMS', 'INT_ARRAY') # List of custom radii - # Valid for: GePol in MODE=ATOMS + # Valid for: GePol and TsLess in MODE=ATOMS # Valid values: array of doubles cavity.add_kw('RADII', 'DBL_ARRAY') # List of spheres - # Valid for: GePol in MODE=EXPLICIT + # Valid for: GePol and TsLess in MODE=EXPLICIT # Valid values: array of doubles in format [x, y, z, R] cavity.add_kw('SPHERES', 'DBL_ARRAY', callback=verify_spheres) top.add_sect(cavity) @@ -373,7 +383,7 @@ def verify_top(section): def verify_cavity(section): - allowed = ('GEPOL', 'RESTART') + allowed = ('GEPOL', 'TSLESS', 'RESTART') type = section.get('TYPE') if (type.get() not in allowed): print(('Requested {} cavity is not among the allowed types: {}'.format(type, allowed))) @@ -385,6 +395,8 @@ def verify_cavity(section): convert_area_scalar(section['AREA']) if (section['MINRADIUS'].is_set()): convert_length_scalar(section['MINRADIUS']) + if (section['MINDISTANCE'].is_set()): + convert_length_scalar(section['MINDISTANCE']) if (type.get() == 'GEPOL'): area = section.get('AREA') @@ -397,6 +409,12 @@ def verify_cavity(section): if (mr < 0.4): print(('Requested minimal radius for added spheres too small: {}. Minimal value is: 0.4 au'.format(mr))) sys.exit(1) + elif (type.get() == 'TSLESS'): + area = section.get('AREA') + a = area.get() + if (a < 0.01): + print(('Requested area value too small: {}. Minimal value is: 0.01 au^2'.format(a))) + sys.exit(1) elif (type.get() == 'RESTART'): npzfile = section.get('NPZFILE') # Restart string is the filename, with extension, where the cavity specs were saved.