From cee8338effd50f104d57f088c5382b5cab794a4e Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Mon, 9 Oct 2023 14:02:45 +0100 Subject: [PATCH 01/15] WIP: Add Schwenke s gas phase H2O potential --- h2opes.f | 493 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 493 insertions(+) create mode 100644 h2opes.f diff --git a/h2opes.f b/h2opes.f new file mode 100644 index 00000000..57a4df53 --- /dev/null +++ b/h2opes.f @@ -0,0 +1,493 @@ + program h2pes + implicit real*8 (a-h,o-z) + dimension rij(1,3),v(1) + rij(1,1) = 1.8 + rij(1,2) = 1.8 + rij(1,3) = 1.7 + call vibpot(rij, v, 1) + print *, "V = ", v(1) + contains + subroutine vibpot(rij,v,n) +c +c pes for h2o, +c Harry Partridge and David W. Schwenke, J. Chem. Phys., +c submitted Aug. 28, 1996. +c rij(i,1)& rij(i,2) are oh distances in au +c rij(i,3) is hoh angle in rad +c v(i) is pes in au +c n is number of geometries +c + dimension rij(n,3),v(n),c5z(245),cbasis(245),ccore(245), + $ crest(245),idx(245,3),fmat(15,3) + data (idx(i,1),i=1,245)/ + $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, + $ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + $ 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + $ 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, + $ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + $ 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, + $ 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, + $ 6, 6, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, + $ 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 5, 5, + $ 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, + $ 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, + $ 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 9, 9, + $ 9, 9, 9, 9, 9/ + data (idx(i,2),i=1,245)/ + $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + $ 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, + $ 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, + $ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + $ 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, + $ 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, + $ 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, + $ 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, + $ 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, + $ 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, + $ 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 1, 1, + $ 1, 1, 1, 1, 1/ + data (idx(i,3),i=1,245)/ + $ 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15, 1, 2, 3, 4, 5, + $ 6, 7, 8, 9,10,11,12,13,14, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11, + $12,13, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13, 1, 2, 3, 4, 5, + $ 6, 7, 8, 9,10,11,12, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12, 1, + $ 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, + $11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8, + $ 9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 1, 2, 3, 4, 5, 6, 7, 8, + $ 9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, + $ 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, + $ 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, + $ 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, + $ 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, + $ 3, 4, 5, 6, 7/ + data (c5z(i),i=1,245)/ + $ 4.2278462684916D+04, 4.5859382909906D-02, 9.4804986183058D+03, + $ 7.5485566680955D+02, 1.9865052511496D+03, 4.3768071560862D+02, + $ 1.4466054104131D+03, 1.3591924557890D+02,-1.4299027252645D+03, + $ 6.6966329416373D+02, 3.8065088734195D+03,-5.0582552618154D+02, + $-3.2067534385604D+03, 6.9673382568135D+02, 1.6789085874578D+03, + $-3.5387509130093D+03,-1.2902326455736D+04,-6.4271125232353D+03, + $-6.9346876863641D+03,-4.9765266152649D+02,-3.4380943579627D+03, + $ 3.9925274973255D+03,-1.2703668547457D+04,-1.5831591056092D+04, + $ 2.9431777405339D+04, 2.5071411925779D+04,-4.8518811956397D+04, + $-1.4430705306580D+04, 2.5844109323395D+04,-2.3371683301770D+03, + $ 1.2333872678202D+04, 6.6525207018832D+03,-2.0884209672231D+03, + $-6.3008463062877D+03, 4.2548148298119D+04, 2.1561445953347D+04, + $-1.5517277060400D+05, 2.9277086555691D+04, 2.6154026873478D+05, + $-1.3093666159230D+05,-1.6260425387088D+05, 1.2311652217133D+05, + $-5.1764697159603D+04, 2.5287599662992D+03, 3.0114701659513D+04, + $-2.0580084492150D+03, 3.3617940269402D+04, 1.3503379582016D+04, + $-1.0401149481887D+05,-6.3248258344140D+04, 2.4576697811922D+05, + $ 8.9685253338525D+04,-2.3910076031416D+05,-6.5265145723160D+04, + $ 8.9184290973880D+04,-8.0850272976101D+03,-3.1054961140464D+04, + $-1.3684354599285D+04, 9.3754012976495D+03,-7.4676475789329D+04, + $-1.8122270942076D+05, 2.6987309391410D+05, 4.0582251904706D+05, + $-4.7103517814752D+05,-3.6115503974010D+05, 3.2284775325099D+05, + $ 1.3264691929787D+04, 1.8025253924335D+05,-1.2235925565102D+04, + $-9.1363898120735D+03,-4.1294242946858D+04,-3.4995730900098D+04, + $ 3.1769893347165D+05, 2.8395605362570D+05,-1.0784536354219D+06, + $-5.9451106980882D+05, 1.5215430060937D+06, 4.5943167339298D+05, + $-7.9957883936866D+05,-9.2432840622294D+04, 5.5825423140341D+03, + $ 3.0673594098716D+03, 8.7439532014842D+04, 1.9113438435651D+05, + $-3.4306742659939D+05,-3.0711488132651D+05, 6.2118702580693D+05, + $-1.5805976377422D+04,-4.2038045404190D+05, 3.4847108834282D+05, + $-1.3486811106770D+04, 3.1256632170871D+04, 5.3344700235019D+03, + $ 2.6384242145376D+04, 1.2917121516510D+05,-1.3160848301195D+05, + $-4.5853998051192D+05, 3.5760105069089D+05, 6.4570143281747D+05, + $-3.6980075904167D+05,-3.2941029518332D+05,-3.5042507366553D+05, + $ 2.1513919629391D+03, 6.3403845616538D+04, 6.2152822008047D+04, + $-4.8805335375295D+05,-6.3261951398766D+05, 1.8433340786742D+06, + $ 1.4650263449690D+06,-2.9204939728308D+06,-1.1011338105757D+06, + $ 1.7270664922758D+06, 3.4925947462024D+05,-1.9526251371308D+04, + $-3.2271030511683D+04,-3.7601575719875D+05, 1.8295007005531D+05, + $ 1.5005699079799D+06,-1.2350076538617D+06,-1.8221938812193D+06, + $ 1.5438780841786D+06,-3.2729150692367D+03, 1.0546285883943D+04, + $-4.7118461673723D+04,-1.1458551385925D+05, 2.7704588008958D+05, + $ 7.4145816862032D+05,-6.6864945408289D+05,-1.6992324545166D+06, + $ 6.7487333473248D+05, 1.4361670430046D+06,-2.0837555267331D+05, + $ 4.7678355561019D+05,-1.5194821786066D+04,-1.1987249931134D+05, + $ 1.3007675671713D+05, 9.6641544907323D+05,-5.3379849922258D+05, + $-2.4303858824867D+06, 1.5261649025605D+06, 2.0186755858342D+06, + $-1.6429544469130D+06,-1.7921520714752D+04, 1.4125624734639D+04, + $-2.5345006031695D+04, 1.7853375909076D+05,-5.4318156343922D+04, + $-3.6889685715963D+05, 4.2449670705837D+05, 3.5020329799394D+05, + $ 9.3825886484788D+03,-8.0012127425648D+05, 9.8554789856472D+04, + $ 4.9210554266522D+05,-6.4038493953446D+05,-2.8398085766046D+06, + $ 2.1390360019254D+06, 6.3452935017176D+06,-2.3677386290925D+06, + $-3.9697874352050D+06,-1.9490691547041D+04, 4.4213579019433D+04, + $ 1.6113884156437D+05,-7.1247665213713D+05,-1.1808376404616D+06, + $ 3.0815171952564D+06, 1.3519809705593D+06,-3.4457898745450D+06, + $ 2.0705775494050D+05,-4.3778169926622D+05, 8.7041260169714D+03, + $ 1.8982512628535D+05,-2.9708215504578D+05,-8.8213012222074D+05, + $ 8.6031109049755D+05, 1.0968800857081D+06,-1.0114716732602D+06, + $ 1.9367263614108D+05, 2.8678295007137D+05,-9.4347729862989D+04, + $ 4.4154039394108D+04, 5.3686756196439D+05, 1.7254041770855D+05, + $-2.5310674462399D+06,-2.0381171865455D+06, 3.3780796258176D+06, + $ 7.8836220768478D+05,-1.5307728782887D+05,-3.7573362053757D+05, + $ 1.0124501604626D+06, 2.0929686545723D+06,-5.7305706586465D+06, + $-2.6200352535413D+06, 7.1543745536691D+06,-1.9733601879064D+04, + $ 8.5273008477607D+04, 6.1062454495045D+04,-2.2642508675984D+05, + $ 2.4581653864150D+05,-9.0376851105383D+05,-4.4367930945690D+05, + $ 1.5740351463593D+06, 2.4563041445249D+05,-3.4697646046367D+03, + $-2.1391370322552D+05, 4.2358948404842D+05, 5.6270081955003D+05, + $-8.5007851251980D+05,-6.1182429537130D+05, 5.6690751824341D+05, + $-3.5617502919487D+05,-8.1875263381402D+02,-2.4506258140060D+05, + $ 2.5830513731509D+05, 6.0646114465433D+05,-6.9676584616955D+05, + $ 5.1937406389690D+05, 1.7261913546007D+05,-1.7405787307472D+04, + $-3.8301842660567D+05, 5.4227693205154D+05, 2.5442083515211D+06, + $-1.1837755702370D+06,-1.9381959088092D+06,-4.0642141553575D+05, + $ 1.1840693827934D+04,-1.5334500255967D+05, 4.9098619510989D+05, + $ 6.1688992640977D+05, 2.2351144690009D+05,-1.8550462739570D+06, + $ 9.6815110649918D+03,-8.1526584681055D+04,-8.0810433155289D+04, + $ 3.4520506615177D+05, 2.5509863381419D+05,-1.3331224992157D+05, + $-4.3119301071653D+05,-5.9818343115856D+04, 1.7863692414573D+03, + $ 8.9440694919836D+04,-2.5558967650731D+05,-2.2130423988459D+04, + $ 4.4973674518316D+05,-2.2094939343618D+05/ + data (cbasis(i),i=1,245)/ + $ 6.9770019624764D-04,-2.4209870001642D+01, 1.8113927151562D+01, + $ 3.5107416275981D+01,-5.4600021126735D+00,-4.8731149608386D+01, + $ 3.6007189184766D+01, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-7.7178474355102D+01,-3.8460795013977D+01,-4.6622480912340D+01, + $ 5.5684951167513D+01, 1.2274939911242D+02,-1.4325154752086D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00,-6.0800589055949D+00, + $ 8.6171499453475D+01,-8.4066835441327D+01,-5.8228085624620D+01, + $ 2.0237393793875D+02, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 3.3525582670313D+02, 7.0056962392208D+01,-4.5312502936708D+01, + $-3.0441141194247D+02, 2.8111438108965D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-1.2983583774779D+02, 3.9781671212935D+01, + $-6.6793945229609D+01,-1.9259805675433D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-8.2855757669957D+02,-5.7003072730941D+01, + $-3.5604806670066D+01, 9.6277766002709D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 8.8645622149112D+02,-7.6908409772041D+01, + $ 6.8111763314154D+01, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 2.5090493428062D+02,-2.3622141780572D+02, 5.8155647658455D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 2.8919570295095D+03, + $-1.7871014635921D+02,-1.3515667622500D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-3.6965613754734D+03, 2.1148158286617D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00,-1.4795670139431D+03, + $ 3.6210798138768D+02, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-5.3552886800881D+03, 3.1006384016202D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 1.6241824368764D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 4.3764909606382D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 1.0940849243716D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 3.0743267832931D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00/ + data (ccore(i),i=1,245)/ + $ 2.4332191647159D-02,-2.9749090113656D+01, 1.8638980892831D+01, + $-6.1272361746520D+00, 2.1567487597605D+00,-1.5552044084945D+01, + $ 8.9752150543954D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-3.5693557878741D+02,-3.0398393196894D+00,-6.5936553294576D+00, + $ 1.6056619388911D+01, 7.8061422868204D+01,-8.6270891686359D+01, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00,-3.1688002530217D+01, + $ 3.7586725583944D+01,-3.2725765966657D+01,-5.6458213299259D+00, + $ 2.1502613314595D+01, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 5.2789943583277D+02,-4.2461079404962D+00,-2.4937638543122D+01, + $-1.1963809321312D+02, 2.0240663228078D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-6.2574211352272D+02,-6.9617539465382D+00, + $-5.9440243471241D+01, 1.4944220180218D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-1.2851139918332D+03,-6.5043516710835D+00, + $ 4.0410829440249D+01,-6.7162452402027D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 1.0031942127832D+03, 7.6137226541944D+01, + $-2.7279242226902D+01, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-3.3059000871075D+01, 2.4384498749480D+01,-1.4597931874215D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 1.6559579606045D+03, + $ 1.5038996611400D+02,-7.3865347730818D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-1.9738401290808D+03,-1.4149993809415D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00,-1.2756627454888D+02, + $ 4.1487702227579D+01, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-1.7406770966429D+03,-9.3812204399266D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-1.1890301282216D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 2.3723447727360D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-1.0279968223292D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 5.7153838472603D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00/ + data (crest(i),i=1,245)/ + $ 0.0000000000000D+00,-4.7430930170000D+00,-1.4422132560000D+01, + $-1.8061146510000D+01, 7.5186735000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-2.7962099800000D+02, 1.7616414260000D+01,-9.9741392630000D+01, + $ 7.1402447000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00,-7.8571336480000D+01, + $ 5.2434353250000D+01, 7.7696745000000D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 1.7799123760000D+02, 1.4564532380000D+02, 2.2347226000000D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-4.3823284100000D+02,-7.2846553000000D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-2.6752313750000D+02, 3.6170310000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00/ + data reoh,thetae,b1,roh,alphaoh,deoh,phh1,phh2/0.958649d0, + $ 104.3475d0,2.0d0,0.9519607159623009d0,2.587949757553683d0, + $ 42290.92019288289d0,16.94879431193463d0,12.66426998162947d0/ + data f5z,fbasis,fcore,frest/0.99967788500000d0, + $ 0.15860145369897d0,-1.6351695982132d0,1d0/ + save + data ifirst/0/ + if(ifirst.eq.0)then + ifirst=1 + write(6,1) + 1 format(/1x,'pes for h2o', + $ /1x,'by Harry Partridge and David W. Schwenke', + $ /1x,'submitted to J. Chem. Phys. Aug. 28, 1996') + write(6,56) + 56 format(/1x,'parameters before adjustment') + write(6,55)phh1,phh2,deoh,alphaoh,roh + 55 format(/1x,'two body potential parameters:', + $ /1x,'hh: phh1 = ',f10.1,' phh2 = ',f5.2, + $ /1x,'oh: deoh = ',f10.1,' alpha = ',f7.4, + $ ' re = ',f7.4) + write(6,4)reoh,thetae,b1 + 4 format(/1x,'three body parameters:', + $ /1x,'reoh = ',f10.4,' thetae = ',f10.4, + $ /1x,'betaoh = ',f10.4, + $ /1x,' i j k',7x,'c5z',9x,'cbasis',10x,'ccore', + $ 10x,'crest') + do 2 i=1,245 + write(6,5)(idx(i,j)-1,j=1,3),c5z(i),cbasis(i),ccore(i),crest(i) + 5 format(1x,3i5,1p4e15.7) + c5z(i)=f5z*c5z(i)+fbasis*cbasis(i)+fcore*ccore(i) + $ +frest*crest(i) + 2 continue + write(6,57)f5z,fbasis,fcore,frest + 57 format(/1x,'adjusting parameters using scale factors ', + $ /1x,'f5z = ',f11.8, + $ /1x,'fbasis = ',f11.8, + $ /1x,'fcore = ',f11.8, + $ /1x,'frest = ',f11.8) + phh1=phh1*f5z + deoh=deoh*f5z + write(6,55)phh1,phh2,deoh,alphaoh,roh + write(6,58)reoh,thetae,b1,((idx(i,j)-1,j=1,3),c5z(i),i=1,245) + 58 format(/1x,'three body parameters:', + $ /1x,'reoh = ',f10.4,' thetae = ',f10.4, + $ /1x,'betaoh = ',f10.4, + $ /1x,' i j k cijk', + $ /(1x,3i5,1pe15.7)) +c +c convert parameters from 1/cm, angstrom to a.u. +c + reoh=reoh/0.529177249d0 + b1=b1*0.529177249d0*0.529177249d0 + do 3 i=1,245 + c5z(i)=c5z(i)*4.556335d-6 + 3 continue + rad=acos(-1d0)/1.8d2 + ce=cos(thetae*rad) + phh1=phh1*exp(phh2) + phh1=phh1*4.556335d-6 + phh2=phh2*0.529177249d0 + deoh=deoh*4.556335d-6 + roh=roh/0.529177249d0 + alphaoh=alphaoh*0.529177249d0 + c5z(1)=c5z(1)*2d0 + end if + do 6 i=1,n + x1=(rij(i,1)-reoh)/reoh + x2=(rij(i,2)-reoh)/reoh + x3=cos(rij(i,3))-ce + rhh=sqrt(rij(i,1)**2+rij(i,2)**2 + $ -2d0*rij(i,1)*rij(i,2)*cos(rij(i,3))) + vhh=phh1*exp(-phh2*rhh) + ex=exp(-alphaoh*(rij(i,1)-roh)) + voh1=deoh*ex*(ex-2d0) + ex=exp(-alphaoh*(rij(i,2)-roh)) + voh2=deoh*ex*(ex-2d0) + fmat(1,1)=1d0 + fmat(1,2)=1d0 + fmat(1,3)=1d0 + do 10 j=2,15 + fmat(j,1)=fmat(j-1,1)*x1 + fmat(j,2)=fmat(j-1,2)*x2 + fmat(j,3)=fmat(j-1,3)*x3 + 10 continue + v(i)=0d0 + do 12 j=2,245 + term=c5z(j)*(fmat(idx(j,1),1)*fmat(idx(j,2),2) + $ +fmat(idx(j,2),1)*fmat(idx(j,1),2)) + $ *fmat(idx(j,3),3) + v(i)=v(i)+term + 12 continue + v(i)=v(i)*exp(-b1*((rij(i,1)-reoh)**2+(rij(i,2)-reoh)**2)) + $ +c5z(1) + $ +voh1+voh2+vhh + 6 continue + return + end + end From d2130db5ee31212b69f29b2b5c0d3db89e348107 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Thu, 25 Jan 2024 15:34:54 +0000 Subject: [PATCH 02/15] Integrate the Schwenke potential --- .fprettify.rc | 2 +- h2opes.f => h2o_schwenke_orig.f | 0 src/Makefile | 5 +- src/force_abin.F90 | 4 +- src/force_h2o.F90 | 53 ++++ src/forces.F90 | 3 + src/fortran_interfaces.F90 | 8 + src/h2o_schwenke.f | 486 ++++++++++++++++++++++++++++++++ src/utils.F90 | 4 +- 9 files changed, 559 insertions(+), 6 deletions(-) rename h2opes.f => h2o_schwenke_orig.f (100%) create mode 100644 src/force_h2o.F90 create mode 100644 src/h2o_schwenke.f diff --git a/.fprettify.rc b/.fprettify.rc index 84df3017..7ad662da 100644 --- a/.fprettify.rc +++ b/.fprettify.rc @@ -26,4 +26,4 @@ enable-decl=True disable-fypp=True case=[1,1,1,2] -exclude=[random.F90,fftw3.F90,force_cp2k.F90] +exclude=[random.F90,fftw3.F90,force_cp2k.F90, h2o_schwenke.f] diff --git a/h2opes.f b/h2o_schwenke_orig.f similarity index 100% rename from h2opes.f rename to h2o_schwenke_orig.f diff --git a/src/Makefile b/src/Makefile index eb072665..ff25d084 100755 --- a/src/Makefile +++ b/src/Makefile @@ -1,6 +1,6 @@ F_OBJS := constants.o fortran_interfaces.o error.o modules.o mpi_wrapper.o files.o utils.o io.o random.o arrays.o qmmm.o fftw_interface.o \ shake.o nosehoover.o gle.o transform.o potentials.o force_spline.o estimators.o ekin.o vinit.o plumed.o \ - remd.o force_bound.o water.o force_cp2k.o sh_integ.o surfacehop.o landau_zener.o\ + remd.o force_bound.o water.o h2o_schwenke.o force_h2o.o force_cp2k.o sh_integ.o surfacehop.o landau_zener.o\ force_mm.o tera_mpi_api.o force_abin.o force_tcpb.o force_tera.o force_terash.o en_restraint.o analyze_ext_template.o geom_analysis.o analysis.o \ minimizer.o mdstep.o forces.o cmdline.o init.o @@ -27,5 +27,8 @@ clean: %.o: %.F90 $(FC) $(FFLAGS) $(DFLAGS) -c $< +%.o: %.f + $(FC) $(FFLAGS) $(DFLAGS) -c $< + %.o: %.cpp $(CXX) $(CXXFLAGS) $(DFLAGS) -c $< diff --git a/src/force_abin.F90 b/src/force_abin.F90 index 26d03f0c..38b96f95 100644 --- a/src/force_abin.F90 +++ b/src/force_abin.F90 @@ -5,7 +5,7 @@ ! The basic workflow is very simple: ! 1. ABIN writes current geometry into file geom.dat ! in a XYZ format without the header. -! 2. ABIN launches the shell script POT/r.pot +! 2. ABIN launches the shell script /r. ! 3. The shellscript does a few things: ! i) Takes the input geometry and prepares the input file ! ii) Launches the QM program @@ -15,7 +15,7 @@ ! ! NOTE: We append bead index to every file name so that we can ! call the interface in parallel in PIMD simulations. -! NOTE: Interface for Surface Hopping is a bit more complicated, +! NOTE: The interface for Surface Hopping is a bit more complicated, ! see interfaces/MOLPRO-SH/r.molpro-sh module mod_shell_interface_private use mod_const, only: DP, ANG diff --git a/src/force_h2o.F90 b/src/force_h2o.F90 new file mode 100644 index 00000000..b33b61b9 --- /dev/null +++ b/src/force_h2o.F90 @@ -0,0 +1,53 @@ +! This is invoked via pot='_h2o_' +module mod_h2o_pes + use mod_const, only: DP + use mod_error, only: fatal_error + implicit none + private + public :: force_h2o_schwenke + save +contains + + subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, walkmax) + use mod_utils, only: get_distance, get_angle + use mod_const, only: PI + real(DP), intent(in) :: x(:, :), y(:, :), z(:, :) + real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :) + real(DP), intent(inout) :: eclas + integer, intent(in) :: natom, walkmax + ! Internal water coordinates + real(DP) :: rOH1, rOH2, aHOH_deg, aHOH_rad + real(DP) :: rij(walkmax, 3) + real(DP) :: Epot(walkmax) + integer :: iw + + ! TODO: Check atom names + if (natom /= 3) then + call fatal_error(__FILE__, __LINE__, "This is not a water molecule!") + end if + + ! Compute internal coordinates from cartesians + ! OH distances in bohrs, HOH angle in radians + do iw = 1, walkmax + rOH1 = get_distance(x, y, z, 1, 2, iw) + rOH2 = get_distance(x, y, z, 1, 3, iw) + aHOH_deg = get_angle(x, y, z, 2, 1, 3, iw) + aHOH_rad = aHOH_deg * PI / 180.0D0 + print*, 'Angle = ', aHOH_deg, aHOH_rad + rij(iw, 1) = rOH1 + rij(iw, 2) = rOH2 + rij(iw, 3) = aHOH_rad + end do + + call pes_h2o_schwenke(rij, Epot, walkmax) + + ! TODO: Implement numerical forces + fx = 0.0D0; fy = 0.0D0; fz = 0.0D0 + + do iw = 1, walkmax + eclas = eclas + Epot(iw) + end do + eclas = eclas / walkmax + end subroutine force_h2o_schwenke + +end module mod_h2o_pes diff --git a/src/forces.F90 b/src/forces.F90 index 3eeb9abb..7a7478d3 100644 --- a/src/forces.F90 +++ b/src/forces.F90 @@ -128,6 +128,7 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, e_pot, chpot, walkmax) use mod_const, only: DP use mod_general, only: natom, ipimd use mod_water, only: watpot + use mod_h2o_pes, only: force_h2o_schwenke use mod_force_mm, only: force_mm use mod_sbc, only: force_sbc, isbc use mod_plumed, only: iplumed, force_plumed @@ -155,6 +156,8 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, e_pot, chpot, walkmax) call force_mm(x, y, z, fx, fy, fz, eclas, walkmax) case ("_mmwater_") call force_water(x, y, z, fx, fy, fz, eclas, natom, walkmax, watpot) + case ("_h2o_") + call force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, walkmax) case ("_splined_grid_") ! Only 1D spline grid supported at the moment call force_splined_grid(x, fx, eclas, walkmax) diff --git a/src/fortran_interfaces.F90 b/src/fortran_interfaces.F90 index 8352ab32..2d754359 100644 --- a/src/fortran_interfaces.F90 +++ b/src/fortran_interfaces.F90 @@ -63,6 +63,14 @@ function usleep(useconds) bind(c, name='usleep') integer(kind=C_INT) :: usleep end function usleep + ! Returns a potential energy of a water molecule + ! using Schwenke potential, see h2o_schwenke.f + subroutine pes_h2o_schwenke(rij, v, n) + import :: DP + integer, intent(in) :: n + real(DP) :: rij(n, 3), v(n) + end subroutine pes_h2o_schwenke + end interface end module mod_interfaces diff --git a/src/h2o_schwenke.f b/src/h2o_schwenke.f new file mode 100644 index 00000000..4a5b15de --- /dev/null +++ b/src/h2o_schwenke.f @@ -0,0 +1,486 @@ + subroutine pes_h2o_schwenke(rij, v, n) + implicit real*8 (a-h,o-z) + integer, intent(in) :: n + integer :: i, idx, j, ifirst +c +c pes for h2o, +c Harry Partridge and David W. Schwenke, J. Chem. Phys., +c submitted Aug. 28, 1996. +c rij(i,1)& rij(i,2) are oh distances in au +c rij(i,3) is hoh angle in rad +c v(i) is pes in au +c n is number of geometries +c + dimension rij(n,3),v(n),c5z(245),cbasis(245),ccore(245), + $ crest(245),idx(245,3),fmat(15,3) + data (idx(i,1),i=1,245)/ + $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, + $ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + $ 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + $ 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, + $ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + $ 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, + $ 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, + $ 6, 6, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, + $ 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 5, 5, + $ 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, + $ 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, + $ 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 9, 9, + $ 9, 9, 9, 9, 9/ + data (idx(i,2),i=1,245)/ + $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + $ 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, + $ 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, + $ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + $ 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, + $ 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, + $ 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, + $ 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, + $ 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, + $ 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, + $ 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 1, 1, + $ 1, 1, 1, 1, 1/ + data (idx(i,3),i=1,245)/ + $ 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15, 1, 2, 3, 4, 5, + $ 6, 7, 8, 9,10,11,12,13,14, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11, + $12,13, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13, 1, 2, 3, 4, 5, + $ 6, 7, 8, 9,10,11,12, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12, 1, + $ 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, + $11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8, + $ 9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 1, 2, 3, 4, 5, 6, 7, 8, + $ 9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, + $ 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, + $ 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, + $ 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, + $ 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, + $ 3, 4, 5, 6, 7/ + data (c5z(i),i=1,245)/ + $ 4.2278462684916D+04, 4.5859382909906D-02, 9.4804986183058D+03, + $ 7.5485566680955D+02, 1.9865052511496D+03, 4.3768071560862D+02, + $ 1.4466054104131D+03, 1.3591924557890D+02,-1.4299027252645D+03, + $ 6.6966329416373D+02, 3.8065088734195D+03,-5.0582552618154D+02, + $-3.2067534385604D+03, 6.9673382568135D+02, 1.6789085874578D+03, + $-3.5387509130093D+03,-1.2902326455736D+04,-6.4271125232353D+03, + $-6.9346876863641D+03,-4.9765266152649D+02,-3.4380943579627D+03, + $ 3.9925274973255D+03,-1.2703668547457D+04,-1.5831591056092D+04, + $ 2.9431777405339D+04, 2.5071411925779D+04,-4.8518811956397D+04, + $-1.4430705306580D+04, 2.5844109323395D+04,-2.3371683301770D+03, + $ 1.2333872678202D+04, 6.6525207018832D+03,-2.0884209672231D+03, + $-6.3008463062877D+03, 4.2548148298119D+04, 2.1561445953347D+04, + $-1.5517277060400D+05, 2.9277086555691D+04, 2.6154026873478D+05, + $-1.3093666159230D+05,-1.6260425387088D+05, 1.2311652217133D+05, + $-5.1764697159603D+04, 2.5287599662992D+03, 3.0114701659513D+04, + $-2.0580084492150D+03, 3.3617940269402D+04, 1.3503379582016D+04, + $-1.0401149481887D+05,-6.3248258344140D+04, 2.4576697811922D+05, + $ 8.9685253338525D+04,-2.3910076031416D+05,-6.5265145723160D+04, + $ 8.9184290973880D+04,-8.0850272976101D+03,-3.1054961140464D+04, + $-1.3684354599285D+04, 9.3754012976495D+03,-7.4676475789329D+04, + $-1.8122270942076D+05, 2.6987309391410D+05, 4.0582251904706D+05, + $-4.7103517814752D+05,-3.6115503974010D+05, 3.2284775325099D+05, + $ 1.3264691929787D+04, 1.8025253924335D+05,-1.2235925565102D+04, + $-9.1363898120735D+03,-4.1294242946858D+04,-3.4995730900098D+04, + $ 3.1769893347165D+05, 2.8395605362570D+05,-1.0784536354219D+06, + $-5.9451106980882D+05, 1.5215430060937D+06, 4.5943167339298D+05, + $-7.9957883936866D+05,-9.2432840622294D+04, 5.5825423140341D+03, + $ 3.0673594098716D+03, 8.7439532014842D+04, 1.9113438435651D+05, + $-3.4306742659939D+05,-3.0711488132651D+05, 6.2118702580693D+05, + $-1.5805976377422D+04,-4.2038045404190D+05, 3.4847108834282D+05, + $-1.3486811106770D+04, 3.1256632170871D+04, 5.3344700235019D+03, + $ 2.6384242145376D+04, 1.2917121516510D+05,-1.3160848301195D+05, + $-4.5853998051192D+05, 3.5760105069089D+05, 6.4570143281747D+05, + $-3.6980075904167D+05,-3.2941029518332D+05,-3.5042507366553D+05, + $ 2.1513919629391D+03, 6.3403845616538D+04, 6.2152822008047D+04, + $-4.8805335375295D+05,-6.3261951398766D+05, 1.8433340786742D+06, + $ 1.4650263449690D+06,-2.9204939728308D+06,-1.1011338105757D+06, + $ 1.7270664922758D+06, 3.4925947462024D+05,-1.9526251371308D+04, + $-3.2271030511683D+04,-3.7601575719875D+05, 1.8295007005531D+05, + $ 1.5005699079799D+06,-1.2350076538617D+06,-1.8221938812193D+06, + $ 1.5438780841786D+06,-3.2729150692367D+03, 1.0546285883943D+04, + $-4.7118461673723D+04,-1.1458551385925D+05, 2.7704588008958D+05, + $ 7.4145816862032D+05,-6.6864945408289D+05,-1.6992324545166D+06, + $ 6.7487333473248D+05, 1.4361670430046D+06,-2.0837555267331D+05, + $ 4.7678355561019D+05,-1.5194821786066D+04,-1.1987249931134D+05, + $ 1.3007675671713D+05, 9.6641544907323D+05,-5.3379849922258D+05, + $-2.4303858824867D+06, 1.5261649025605D+06, 2.0186755858342D+06, + $-1.6429544469130D+06,-1.7921520714752D+04, 1.4125624734639D+04, + $-2.5345006031695D+04, 1.7853375909076D+05,-5.4318156343922D+04, + $-3.6889685715963D+05, 4.2449670705837D+05, 3.5020329799394D+05, + $ 9.3825886484788D+03,-8.0012127425648D+05, 9.8554789856472D+04, + $ 4.9210554266522D+05,-6.4038493953446D+05,-2.8398085766046D+06, + $ 2.1390360019254D+06, 6.3452935017176D+06,-2.3677386290925D+06, + $-3.9697874352050D+06,-1.9490691547041D+04, 4.4213579019433D+04, + $ 1.6113884156437D+05,-7.1247665213713D+05,-1.1808376404616D+06, + $ 3.0815171952564D+06, 1.3519809705593D+06,-3.4457898745450D+06, + $ 2.0705775494050D+05,-4.3778169926622D+05, 8.7041260169714D+03, + $ 1.8982512628535D+05,-2.9708215504578D+05,-8.8213012222074D+05, + $ 8.6031109049755D+05, 1.0968800857081D+06,-1.0114716732602D+06, + $ 1.9367263614108D+05, 2.8678295007137D+05,-9.4347729862989D+04, + $ 4.4154039394108D+04, 5.3686756196439D+05, 1.7254041770855D+05, + $-2.5310674462399D+06,-2.0381171865455D+06, 3.3780796258176D+06, + $ 7.8836220768478D+05,-1.5307728782887D+05,-3.7573362053757D+05, + $ 1.0124501604626D+06, 2.0929686545723D+06,-5.7305706586465D+06, + $-2.6200352535413D+06, 7.1543745536691D+06,-1.9733601879064D+04, + $ 8.5273008477607D+04, 6.1062454495045D+04,-2.2642508675984D+05, + $ 2.4581653864150D+05,-9.0376851105383D+05,-4.4367930945690D+05, + $ 1.5740351463593D+06, 2.4563041445249D+05,-3.4697646046367D+03, + $-2.1391370322552D+05, 4.2358948404842D+05, 5.6270081955003D+05, + $-8.5007851251980D+05,-6.1182429537130D+05, 5.6690751824341D+05, + $-3.5617502919487D+05,-8.1875263381402D+02,-2.4506258140060D+05, + $ 2.5830513731509D+05, 6.0646114465433D+05,-6.9676584616955D+05, + $ 5.1937406389690D+05, 1.7261913546007D+05,-1.7405787307472D+04, + $-3.8301842660567D+05, 5.4227693205154D+05, 2.5442083515211D+06, + $-1.1837755702370D+06,-1.9381959088092D+06,-4.0642141553575D+05, + $ 1.1840693827934D+04,-1.5334500255967D+05, 4.9098619510989D+05, + $ 6.1688992640977D+05, 2.2351144690009D+05,-1.8550462739570D+06, + $ 9.6815110649918D+03,-8.1526584681055D+04,-8.0810433155289D+04, + $ 3.4520506615177D+05, 2.5509863381419D+05,-1.3331224992157D+05, + $-4.3119301071653D+05,-5.9818343115856D+04, 1.7863692414573D+03, + $ 8.9440694919836D+04,-2.5558967650731D+05,-2.2130423988459D+04, + $ 4.4973674518316D+05,-2.2094939343618D+05/ + data (cbasis(i),i=1,245)/ + $ 6.9770019624764D-04,-2.4209870001642D+01, 1.8113927151562D+01, + $ 3.5107416275981D+01,-5.4600021126735D+00,-4.8731149608386D+01, + $ 3.6007189184766D+01, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-7.7178474355102D+01,-3.8460795013977D+01,-4.6622480912340D+01, + $ 5.5684951167513D+01, 1.2274939911242D+02,-1.4325154752086D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00,-6.0800589055949D+00, + $ 8.6171499453475D+01,-8.4066835441327D+01,-5.8228085624620D+01, + $ 2.0237393793875D+02, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 3.3525582670313D+02, 7.0056962392208D+01,-4.5312502936708D+01, + $-3.0441141194247D+02, 2.8111438108965D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-1.2983583774779D+02, 3.9781671212935D+01, + $-6.6793945229609D+01,-1.9259805675433D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-8.2855757669957D+02,-5.7003072730941D+01, + $-3.5604806670066D+01, 9.6277766002709D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 8.8645622149112D+02,-7.6908409772041D+01, + $ 6.8111763314154D+01, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 2.5090493428062D+02,-2.3622141780572D+02, 5.8155647658455D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 2.8919570295095D+03, + $-1.7871014635921D+02,-1.3515667622500D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-3.6965613754734D+03, 2.1148158286617D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00,-1.4795670139431D+03, + $ 3.6210798138768D+02, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-5.3552886800881D+03, 3.1006384016202D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 1.6241824368764D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 4.3764909606382D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 1.0940849243716D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 3.0743267832931D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00/ + data (ccore(i),i=1,245)/ + $ 2.4332191647159D-02,-2.9749090113656D+01, 1.8638980892831D+01, + $-6.1272361746520D+00, 2.1567487597605D+00,-1.5552044084945D+01, + $ 8.9752150543954D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-3.5693557878741D+02,-3.0398393196894D+00,-6.5936553294576D+00, + $ 1.6056619388911D+01, 7.8061422868204D+01,-8.6270891686359D+01, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00,-3.1688002530217D+01, + $ 3.7586725583944D+01,-3.2725765966657D+01,-5.6458213299259D+00, + $ 2.1502613314595D+01, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 5.2789943583277D+02,-4.2461079404962D+00,-2.4937638543122D+01, + $-1.1963809321312D+02, 2.0240663228078D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-6.2574211352272D+02,-6.9617539465382D+00, + $-5.9440243471241D+01, 1.4944220180218D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-1.2851139918332D+03,-6.5043516710835D+00, + $ 4.0410829440249D+01,-6.7162452402027D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 1.0031942127832D+03, 7.6137226541944D+01, + $-2.7279242226902D+01, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-3.3059000871075D+01, 2.4384498749480D+01,-1.4597931874215D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 1.6559579606045D+03, + $ 1.5038996611400D+02,-7.3865347730818D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-1.9738401290808D+03,-1.4149993809415D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00,-1.2756627454888D+02, + $ 4.1487702227579D+01, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-1.7406770966429D+03,-9.3812204399266D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-1.1890301282216D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 2.3723447727360D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-1.0279968223292D+03, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 5.7153838472603D+02, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00/ + data (crest(i),i=1,245)/ + $ 0.0000000000000D+00,-4.7430930170000D+00,-1.4422132560000D+01, + $-1.8061146510000D+01, 7.5186735000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $-2.7962099800000D+02, 1.7616414260000D+01,-9.9741392630000D+01, + $ 7.1402447000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00,-7.8571336480000D+01, + $ 5.2434353250000D+01, 7.7696745000000D+01, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 1.7799123760000D+02, 1.4564532380000D+02, 2.2347226000000D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-4.3823284100000D+02,-7.2846553000000D+02, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00,-2.6752313750000D+02, 3.6170310000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00/ + data reoh,thetae,b1,roh,alphaoh,deoh,phh1,phh2/0.958649d0, + $ 104.3475d0,2.0d0,0.9519607159623009d0,2.587949757553683d0, + $ 42290.92019288289d0,16.94879431193463d0,12.66426998162947d0/ + data f5z,fbasis,fcore,frest/0.99967788500000d0, + $ 0.15860145369897d0,-1.6351695982132d0,1d0/ + save + data ifirst/0/ + if(ifirst.eq.0)then + ifirst=1 + write(6,1) + 1 format(/1x,'pes for h2o', + $ /1x,'by Harry Partridge and David W. Schwenke', + $ /1x,'submitted to J. Chem. Phys. Aug. 28, 1996') + write(6,56) + 56 format(/1x,'parameters before adjustment') + write(6,55)phh1,phh2,deoh,alphaoh,roh + 55 format(/1x,'two body potential parameters:', + $ /1x,'hh: phh1 = ',f10.1,' phh2 = ',f5.2, + $ /1x,'oh: deoh = ',f10.1,' alpha = ',f7.4, + $ ' re = ',f7.4) + write(6,4)reoh,thetae,b1 + 4 format(/1x,'three body parameters:', + $ /1x,'reoh = ',f10.4,' thetae = ',f10.4, + $ /1x,'betaoh = ',f10.4, + $ /1x,' i j k',7x,'c5z',9x,'cbasis',10x,'ccore', + $ 10x,'crest') + do 2 i=1,245 +c write(6,5)(idx(i,j)-1,j=1,3),c5z(i),cbasis(i),ccore(i),crest(i) +c 5 format(1x,3i5,1p4e15.7) + c5z(i)=f5z*c5z(i)+fbasis*cbasis(i)+fcore*ccore(i) + $ +frest*crest(i) + 2 continue + write(6,57)f5z,fbasis,fcore,frest + 57 format(/1x,'adjusting parameters using scale factors ', + $ /1x,'f5z = ',f11.8, + $ /1x,'fbasis = ',f11.8, + $ /1x,'fcore = ',f11.8, + $ /1x,'frest = ',f11.8) + phh1=phh1*f5z + deoh=deoh*f5z + write(6,55)phh1,phh2,deoh,alphaoh,roh +c write(6,58)reoh,thetae,b1,((idx(i,j)-1,j=1,3),c5z(i),i=1,245) +c 58 format(/1x,'three body parameters:', +c $ /1x,'reoh = ',f10.4,' thetae = ',f10.4, +c $ /1x,'betaoh = ',f10.4, +c $ /1x,' i j k cijk', +c $ /(1x,3i5,1pe15.7)) +c +c convert parameters from 1/cm, angstrom to a.u. +c + reoh=reoh/0.529177249d0 + b1=b1*0.529177249d0*0.529177249d0 + do 3 i=1,245 + c5z(i)=c5z(i)*4.556335d-6 + 3 continue + rad=acos(-1d0)/1.8d2 + ce=cos(thetae*rad) + phh1=phh1*exp(phh2) + phh1=phh1*4.556335d-6 + phh2=phh2*0.529177249d0 + deoh=deoh*4.556335d-6 + roh=roh/0.529177249d0 + alphaoh=alphaoh*0.529177249d0 + c5z(1)=c5z(1)*2d0 + end if + do 6 i=1,n + x1=(rij(i,1)-reoh)/reoh + x2=(rij(i,2)-reoh)/reoh + x3=cos(rij(i,3))-ce + rhh=sqrt(rij(i,1)**2+rij(i,2)**2 + $ -2d0*rij(i,1)*rij(i,2)*cos(rij(i,3))) + vhh=phh1*exp(-phh2*rhh) + ex=exp(-alphaoh*(rij(i,1)-roh)) + voh1=deoh*ex*(ex-2d0) + ex=exp(-alphaoh*(rij(i,2)-roh)) + voh2=deoh*ex*(ex-2d0) + fmat(1,1)=1d0 + fmat(1,2)=1d0 + fmat(1,3)=1d0 + do 10 j=2,15 + fmat(j,1)=fmat(j-1,1)*x1 + fmat(j,2)=fmat(j-1,2)*x2 + fmat(j,3)=fmat(j-1,3)*x3 + 10 continue + v(i)=0d0 + do 12 j=2,245 + term=c5z(j)*(fmat(idx(j,1),1)*fmat(idx(j,2),2) + $ +fmat(idx(j,2),1)*fmat(idx(j,1),2)) + $ *fmat(idx(j,3),3) + v(i)=v(i)+term + 12 continue + v(i)=v(i)*exp(-b1*((rij(i,1)-reoh)**2+(rij(i,2)-reoh)**2)) + $ +c5z(1) + $ +voh1+voh2+vhh + 6 continue + return + end subroutine diff --git a/src/utils.F90 b/src/utils.F90 index 2dfff0fc..1336be33 100644 --- a/src/utils.F90 +++ b/src/utils.F90 @@ -46,7 +46,7 @@ real(DP) function get_angle(x, y, z, at1, at2, at3, iw) result(angle) vec2x = x(at3, iw) - x(at2, iw) vec2y = y(at3, iw) - y(at2, iw) vec2z = z(at3, iw) - z(at2, iw) - angle = 180 / pi * acos((vec1x * vec2x + vec1y * vec2y + vec1z * vec2z) / & + angle = 180.0D0 / PI * acos((vec1x * vec2x + vec1y * vec2y + vec1z * vec2z) / & & (dsqrt(vec1x**2 + vec1y**2 + vec1z**2) * dsqrt(vec2x**2 + vec2y**2 + vec2z**2))) end function get_angle @@ -81,7 +81,7 @@ real(DP) function get_dihedral(x, y, z, at1, at2, at3, at4, iw, shiftdih) ! TODO: Refactor, make more intermediate results, e.g. ! norms of the normal vectors. ! TODO: Add error handling for malformed dihedral angles to prevent division by zero - get_dihedral = 180 / pi * acos( & + get_dihedral = 180.0D0 / PI * acos( & & (norm1x * norm2x + norm1y * norm2y + norm1z * norm2z) / & & (dsqrt(norm1x**2 + norm1y**2 + norm1z**2) * dsqrt(norm2x**2 + norm2y**2 + norm2z**2)) & & ) From 517f8ed92942d5ec696a2d04946601105b7b11c8 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Thu, 25 Jan 2024 16:35:57 +0000 Subject: [PATCH 03/15] Remove print statement --- src/force_h2o.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/force_h2o.F90 b/src/force_h2o.F90 index b33b61b9..f8bb651a 100644 --- a/src/force_h2o.F90 +++ b/src/force_h2o.F90 @@ -33,7 +33,6 @@ subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, walkmax) rOH2 = get_distance(x, y, z, 1, 3, iw) aHOH_deg = get_angle(x, y, z, 2, 1, 3, iw) aHOH_rad = aHOH_deg * PI / 180.0D0 - print*, 'Angle = ', aHOH_deg, aHOH_rad rij(iw, 1) = rOH1 rij(iw, 2) = rOH2 rij(iw, 3) = aHOH_rad From 4b88ef1a9d01574d2883d906bff2b4bc356a1978 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Thu, 25 Jan 2024 16:36:32 +0000 Subject: [PATCH 04/15] GHA: Setup python-3.11 for autoformat --- .github/workflows/gfortran.yml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/.github/workflows/gfortran.yml b/.github/workflows/gfortran.yml index f91d437f..749697e9 100644 --- a/.github/workflows/gfortran.yml +++ b/.github/workflows/gfortran.yml @@ -39,6 +39,12 @@ jobs: with: fetch-depth: 2 + - name: Set up Python with cache + uses: actions/setup-python@v4 + with: + python-version: '3.11' + cache: pip + - name: Install fprettify dependencies run: pip install --user configargparse From cef431f9c1dbdb2c465f1b2d683c92a0a4e3ef36 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Thu, 25 Jan 2024 16:57:49 +0000 Subject: [PATCH 05/15] No cache --- .github/workflows/gfortran.yml | 7 +------ autoformat.py | 4 ++-- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/.github/workflows/gfortran.yml b/.github/workflows/gfortran.yml index 749697e9..6e2a2e9b 100644 --- a/.github/workflows/gfortran.yml +++ b/.github/workflows/gfortran.yml @@ -36,17 +36,12 @@ jobs: steps: - uses: actions/checkout@v4 - with: - fetch-depth: 2 - - - name: Set up Python with cache uses: actions/setup-python@v4 with: python-version: '3.11' - cache: pip - name: Install fprettify dependencies - run: pip install --user configargparse + run: pip install configargparse - name: Run fprettify run: | diff --git a/autoformat.py b/autoformat.py index 5b3e208e..e42a386d 100755 --- a/autoformat.py +++ b/autoformat.py @@ -31,8 +31,8 @@ try: import configargparse except ImportError: - print("Could not find configargparse package") - print("Please install it e.g. via 'pip install --user configargparse") + print("ERROR: Could not find configargparse package") + print("Please install the package, e.g. 'pip install --user configargparse'") sys.exit(1) if len(sys.argv) == 1: From 2d3d0a7968818f9998326f0d33700ce7995c7126 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Thu, 25 Jan 2024 17:00:59 +0000 Subject: [PATCH 06/15] Use default --- .github/workflows/gfortran.yml | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/.github/workflows/gfortran.yml b/.github/workflows/gfortran.yml index 6e2a2e9b..e2b7d81b 100644 --- a/.github/workflows/gfortran.yml +++ b/.github/workflows/gfortran.yml @@ -29,16 +29,12 @@ env: jobs: format: - name: Format check - runs-on: ubuntu-20.04 - strategy: - fail-fast: false + name: Code format check + runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 uses: actions/setup-python@v4 - with: - python-version: '3.11' - name: Install fprettify dependencies run: pip install configargparse @@ -47,6 +43,7 @@ jobs: run: | ./autoformat.py f=`git ls-files -m`; if [[ -n $f ]];then echo -e "ERROR: Detected unformatted files:\n$f";exit 1;fi + echo "All Fortran files in src/ are formatted correctly!" basic_build: name: Basic build From 78176e7d6d8bb5a317e7bab09ba935b77a48e0af Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Thu, 25 Jan 2024 17:39:36 +0000 Subject: [PATCH 07/15] Whoops --- .github/workflows/gfortran.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/gfortran.yml b/.github/workflows/gfortran.yml index e2b7d81b..6cab5f8f 100644 --- a/.github/workflows/gfortran.yml +++ b/.github/workflows/gfortran.yml @@ -34,7 +34,7 @@ jobs: steps: - uses: actions/checkout@v4 - uses: actions/setup-python@v4 + - uses: actions/setup-python@v4 - name: Install fprettify dependencies run: pip install configargparse From 5063ac3af3897109f752205a872bb5629703e9b4 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Thu, 25 Jan 2024 19:12:05 +0000 Subject: [PATCH 08/15] what is the default python version? --- .github/workflows/gfortran.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/gfortran.yml b/.github/workflows/gfortran.yml index 6cab5f8f..6a146aef 100644 --- a/.github/workflows/gfortran.yml +++ b/.github/workflows/gfortran.yml @@ -34,13 +34,14 @@ jobs: steps: - uses: actions/checkout@v4 - - uses: actions/setup-python@v4 + - uses: actions/setup-python@v5 - name: Install fprettify dependencies run: pip install configargparse - name: Run fprettify run: | + which python && python --version ./autoformat.py f=`git ls-files -m`; if [[ -n $f ]];then echo -e "ERROR: Detected unformatted files:\n$f";exit 1;fi echo "All Fortran files in src/ are formatted correctly!" From e682532aa8e3c084539d712dd0fc43941e7462f9 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Thu, 25 Jan 2024 19:20:40 +0000 Subject: [PATCH 09/15] Py3.11 --- .github/workflows/gfortran.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/gfortran.yml b/.github/workflows/gfortran.yml index 6a146aef..2fc99c90 100644 --- a/.github/workflows/gfortran.yml +++ b/.github/workflows/gfortran.yml @@ -35,13 +35,14 @@ jobs: steps: - uses: actions/checkout@v4 - uses: actions/setup-python@v5 + with: + python-version: '3.11' - name: Install fprettify dependencies run: pip install configargparse - name: Run fprettify run: | - which python && python --version ./autoformat.py f=`git ls-files -m`; if [[ -n $f ]];then echo -e "ERROR: Detected unformatted files:\n$f";exit 1;fi echo "All Fortran files in src/ are formatted correctly!" From 64147ea8aab24339e5e8fe82fae3dcb5ca36752a Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Fri, 26 Jan 2024 18:40:38 +0000 Subject: [PATCH 10/15] Generalize and refactor --- src/force_h2o.F90 | 116 ++++++++++++++++++++++++++++--------- src/forces.F90 | 4 +- src/fortran_interfaces.F90 | 4 +- src/h2o_schwenke.f | 2 +- src/init.F90 | 6 +- 5 files changed, 99 insertions(+), 33 deletions(-) diff --git a/src/force_h2o.F90 b/src/force_h2o.F90 index f8bb651a..e35fa66c 100644 --- a/src/force_h2o.F90 +++ b/src/force_h2o.F90 @@ -1,52 +1,114 @@ -! This is invoked via pot='_h2o_' -module mod_h2o_pes +! Interface to analytical H2O (single water molecule) potentials. +! This is invoked via pot='_h2o_', and the potential is selected via +! h2opot='schwenke' +module mod_force_h2o use mod_const, only: DP + use mod_files, only: stdout, stderr use mod_error, only: fatal_error + use mod_interfaces, only: h2o_pot_schwenke implicit none private - public :: force_h2o_schwenke + public :: initialize_h2o_pot, force_h2o, h2opot + character(len=20) :: h2opot = 'schwenke' save contains - subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, walkmax) - use mod_utils, only: get_distance, get_angle + ! Perform some basic validation + subroutine initialize_h2o_pot(natom, atom_names) + integer, intent(in) :: natom + character(len=2), intent(in) :: atom_names(natom) + + if (natom /= 3) then + call fatal_error(__FILE__, __LINE__, "This is not a water molecule!") + end if + + if (atom_names(1) /= 'O' .or. atom_names(2) /= 'H' .or. atom_names(3) /= 'H') then + write (stderr, *) 'ERROR: Bad element type.' + write (stderr, *) 'Water atoms must be ordered as "O H H"' + call fatal_error(__FILE__, __LINE__, "This is not a water molecule!") + end if + end subroutine initialize_h2o_pot + + ! Compute internal coordinates from cartesians + ! OH distances in bohrs, HOH angle in radians + subroutine get_internal_coords(x, y, z, iw, rOH1, rOH2, aHOH_rad) use mod_const, only: PI + use mod_utils, only: get_distance, get_angle + real(DP), intent(in) :: x(:, :), y(:, :), z(:, :) + integer, intent(in) :: iw + real(DP), intent(out) :: rOH1, rOH2, aHOH_rad + real(DP) :: aHOH_deg + + rOH1 = get_distance(x, y, z, 1, 2, iw) + rOH2 = get_distance(x, y, z, 1, 3, iw) + aHOH_deg = get_angle(x, y, z, 2, 1, 3, iw) + aHOH_rad = aHOH_deg * PI / 180.0D0 + end subroutine get_internal_coords + + ! This is just a wrapper function that selects the H2O potential + ! based on user input. For now only h2opot="schwenke" is supported + subroutine force_h2o(x, y, z, fx, fy, fz, eclas, natom, nbeads) real(DP), intent(in) :: x(:, :), y(:, :), z(:, :) real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :) real(DP), intent(inout) :: eclas - integer, intent(in) :: natom, walkmax - ! Internal water coordinates - real(DP) :: rOH1, rOH2, aHOH_deg, aHOH_rad - real(DP) :: rij(walkmax, 3) - real(DP) :: Epot(walkmax) - integer :: iw + integer, intent(in) :: natom, nbeads - ! TODO: Check atom names - if (natom /= 3) then - call fatal_error(__FILE__, __LINE__, "This is not a water molecule!") + ! TODO: Use function pointers to select the potential and pass it down + if (h2opot == 'schwenke') then + call force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, nbeads) + else + call fatal_error(__FILE__, __LINE__, 'Potential '//trim(h2opot)//' not implemented') end if + end subroutine force_h2o + + subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, nbeads) + real(DP), intent(in) :: x(:, :), y(:, :), z(:, :) + real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :) + real(DP), intent(inout) :: eclas + integer, intent(in) :: natom, nbeads + ! Internal water coordinates + real(DP) :: rOH1, rOH2, aHOH_rad + real(DP) :: rij(nbeads, 3) + real(DP) :: Epot(nbeads) + integer :: iw - ! Compute internal coordinates from cartesians - ! OH distances in bohrs, HOH angle in radians - do iw = 1, walkmax - rOH1 = get_distance(x, y, z, 1, 2, iw) - rOH2 = get_distance(x, y, z, 1, 3, iw) - aHOH_deg = get_angle(x, y, z, 2, 1, 3, iw) - aHOH_rad = aHOH_deg * PI / 180.0D0 + ! The H2O potentials are evaluated in internal coordinates, but ABIN works in cartesians + do iw = 1, nbeads + call get_internal_coords(x, y, z, iw, rOH1, rOH2, aHOH_rad) rij(iw, 1) = rOH1 rij(iw, 2) = rOH2 rij(iw, 3) = aHOH_rad end do - call pes_h2o_schwenke(rij, Epot, walkmax) + ! Calculated the potential for all PI beads at once + call h2o_pot_schwenke(rij, Epot, nbeads) - ! TODO: Implement numerical forces - fx = 0.0D0; fy = 0.0D0; fz = 0.0D0 + call numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads) - do iw = 1, walkmax + ! For Path Integrals, the final energy of the PI necklace + ! is an average over all beads. + ! TODO: Forces need to be appropriately scaled as well! + do iw = 1, nbeads eclas = eclas + Epot(iw) end do - eclas = eclas / walkmax + eclas = eclas / nbeads end subroutine force_h2o_schwenke -end module mod_h2o_pes + ! TODO: Implement numerical forces generally for all potentials + ! For now, they can be implemented here and hardcoded for a specific H2O potential + subroutine numerical_forces(x, y, z, Epot, fx, fy, fz, natom, nbeads) + real(DP), intent(in) :: x(natom, nbeads) + real(DP), intent(in) :: y(natom, nbeads) + real(DP), intent(in) :: z(natom, nbeads) + real(DP), intent(inout) :: fx(natom, nbeads) + real(DP), intent(inout) :: fy(natom, nbeads) + real(DP), intent(inout) :: fz(natom, nbeads) + ! This is the energy for the currrent geometry that has already been calculated + real(DP), intent(in) :: Epot(nbeads) + integer, intent(in) :: natom, nbeads + + ! Need to implement numerical forces first + call fatal_error(__FILE__, __LINE__, 'Numerical forces not yet implemented!') + end subroutine numerical_forces + +end module mod_force_h2o diff --git a/src/forces.F90 b/src/forces.F90 index 7a7478d3..4e527c6e 100644 --- a/src/forces.F90 +++ b/src/forces.F90 @@ -128,7 +128,7 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, e_pot, chpot, walkmax) use mod_const, only: DP use mod_general, only: natom, ipimd use mod_water, only: watpot - use mod_h2o_pes, only: force_h2o_schwenke + use mod_force_h2o, only: force_h2o use mod_force_mm, only: force_mm use mod_sbc, only: force_sbc, isbc use mod_plumed, only: iplumed, force_plumed @@ -157,7 +157,7 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, e_pot, chpot, walkmax) case ("_mmwater_") call force_water(x, y, z, fx, fy, fz, eclas, natom, walkmax, watpot) case ("_h2o_") - call force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, walkmax) + call force_h2o(x, y, z, fx, fy, fz, eclas, natom, walkmax) case ("_splined_grid_") ! Only 1D spline grid supported at the moment call force_splined_grid(x, fx, eclas, walkmax) diff --git a/src/fortran_interfaces.F90 b/src/fortran_interfaces.F90 index 2d754359..696eb911 100644 --- a/src/fortran_interfaces.F90 +++ b/src/fortran_interfaces.F90 @@ -65,11 +65,11 @@ end function usleep ! Returns a potential energy of a water molecule ! using Schwenke potential, see h2o_schwenke.f - subroutine pes_h2o_schwenke(rij, v, n) + subroutine h2o_pot_schwenke(rij, v, n) import :: DP integer, intent(in) :: n real(DP) :: rij(n, 3), v(n) - end subroutine pes_h2o_schwenke + end subroutine h2o_pot_schwenke end interface diff --git a/src/h2o_schwenke.f b/src/h2o_schwenke.f index 4a5b15de..cd63fde9 100644 --- a/src/h2o_schwenke.f +++ b/src/h2o_schwenke.f @@ -1,4 +1,4 @@ - subroutine pes_h2o_schwenke(rij, v, n) + subroutine h2o_pot_schwenke(rij, v, n) implicit real*8 (a-h,o-z) integer, intent(in) :: n integer :: i, idx, j, ifirst diff --git a/src/init.F90 b/src/init.F90 index 14647dd9..52a66e90 100644 --- a/src/init.F90 +++ b/src/init.F90 @@ -45,6 +45,7 @@ subroutine init(dt) use mod_lz use mod_qmmm, only: natqm, natmm use mod_force_mm, only: initialize_mm + use mod_force_h2o, only: initialize_h2o_pot, h2opot use mod_gle use mod_sbc, only: sbc_init, rb_sbc, kb_sbc, isbc, rho use mod_prng_init, only: initialize_prng @@ -126,7 +127,7 @@ subroutine init(dt) namelist /general/ pot, ipimd, mdtype, istage, inormalmodes, nwalk, nstep, icv, ihess, imini, nproc, iqmmm, & nwrite, nwritex, nwritev, nwritef, dt, irandom, nabin, irest, nrest, anal_ext, & isbc, rb_sbc, kb_sbc, gamm, gammthr, conatom, mpi_milisleep, narchive, xyz_units, & - dime, ncalc, idebug, enmini, rho, iknow, watpot, iremd, iplumed, plumedfile, & + dime, ncalc, idebug, enmini, rho, iknow, watpot, h2opot, iremd, iplumed, plumedfile, & en_restraint, en_diff, en_kk, restrain_pot, & pot_ref, nstep_ref, nteraservers, max_mpi_wait_time, cp2k_mpi_beads @@ -511,6 +512,9 @@ subroutine init(dt) if (pot == '_splined_grid_' .or. pot_ref == '_splined_grid_') then call initialize_spline(natom) end if + if (pot == '_h2o_' .or. pot_ref == '_ref_') then + call initialize_h2o_pot(natom, atnames) + end if if (pot == '_mm_' .or. pot_ref == '_mm_') then call initialize_mm(natom, atnames, mm_types, q, LJ_rmin, LJ_eps) end if From 2d67d77d098af9b1789677ce52986f5d5c8fea58 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Fri, 26 Jan 2024 18:50:03 +0000 Subject: [PATCH 11/15] Tiny fix --- .github/workflows/gfortran.yml | 1 - src/init.F90 | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/workflows/gfortran.yml b/.github/workflows/gfortran.yml index 2fc99c90..d7233472 100644 --- a/.github/workflows/gfortran.yml +++ b/.github/workflows/gfortran.yml @@ -45,7 +45,6 @@ jobs: run: | ./autoformat.py f=`git ls-files -m`; if [[ -n $f ]];then echo -e "ERROR: Detected unformatted files:\n$f";exit 1;fi - echo "All Fortran files in src/ are formatted correctly!" basic_build: name: Basic build diff --git a/src/init.F90 b/src/init.F90 index 52a66e90..098b0a68 100644 --- a/src/init.F90 +++ b/src/init.F90 @@ -512,7 +512,7 @@ subroutine init(dt) if (pot == '_splined_grid_' .or. pot_ref == '_splined_grid_') then call initialize_spline(natom) end if - if (pot == '_h2o_' .or. pot_ref == '_ref_') then + if (pot == '_h2o_' .or. pot_ref == '_h2o_') then call initialize_h2o_pot(natom, atnames) end if if (pot == '_mm_' .or. pot_ref == '_mm_') then From 0c39f456d39a797446a1fce1011af0ec08477a4a Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Fri, 26 Jan 2024 18:55:12 +0000 Subject: [PATCH 12/15] Add dummy test --- tests/H2O_SCHWENKE/ERROR.ref | 1 + tests/H2O_SCHWENKE/input.in | 19 +++++++++++++++++++ tests/H2O_SCHWENKE/mini.xyz | 5 +++++ 3 files changed, 25 insertions(+) create mode 100644 tests/H2O_SCHWENKE/ERROR.ref create mode 100644 tests/H2O_SCHWENKE/input.in create mode 100644 tests/H2O_SCHWENKE/mini.xyz diff --git a/tests/H2O_SCHWENKE/ERROR.ref b/tests/H2O_SCHWENKE/ERROR.ref new file mode 100644 index 00000000..783ca168 --- /dev/null +++ b/tests/H2O_SCHWENKE/ERROR.ref @@ -0,0 +1 @@ +ERROR in force_h2o.F90: Numerical forces not yet implemented! diff --git a/tests/H2O_SCHWENKE/input.in b/tests/H2O_SCHWENKE/input.in new file mode 100644 index 00000000..41b29dd9 --- /dev/null +++ b/tests/H2O_SCHWENKE/input.in @@ -0,0 +1,19 @@ +! Test analytical water potential by Schwenke et al + +&general +nstep=1, +irest=0, +idebug=1 + +pot='_h2o_' +h2opot='schwenke' +mdtype='MD', ! classical MD +dt=20., ! number of steps and timestep +nstep=1 +/ + +&nhcopt +inose=0, ! Thermostating: Nose-Hoover 1, microcanonical 0,GLE 2, LE 3 +temp=100 +rem_comrot=.true. ! this is a default value, remove rotations at the beginning +/ diff --git a/tests/H2O_SCHWENKE/mini.xyz b/tests/H2O_SCHWENKE/mini.xyz new file mode 100644 index 00000000..85d61ae3 --- /dev/null +++ b/tests/H2O_SCHWENKE/mini.xyz @@ -0,0 +1,5 @@ +3 + +o 0.000 0.118 0.000 +H 0.757 -0.472 0.000 +h -0.757 -0.472 0.000 From 4b823674e47fa27561febcb541aa10fcbb9244ea Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Fri, 26 Jan 2024 18:56:58 +0000 Subject: [PATCH 13/15] Reorder for test coverage --- src/force_h2o.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/force_h2o.F90 b/src/force_h2o.F90 index e35fa66c..1d7c3667 100644 --- a/src/force_h2o.F90 +++ b/src/force_h2o.F90 @@ -61,7 +61,7 @@ subroutine force_h2o(x, y, z, fx, fy, fz, eclas, natom, nbeads) end if end subroutine force_h2o - subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, nbeads) + subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, Eclas, natom, nbeads) real(DP), intent(in) :: x(:, :), y(:, :), z(:, :) real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :) real(DP), intent(inout) :: eclas @@ -82,16 +82,16 @@ subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, eclas, natom, nbeads) ! Calculated the potential for all PI beads at once call h2o_pot_schwenke(rij, Epot, nbeads) - - call numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads) - ! For Path Integrals, the final energy of the PI necklace ! is an average over all beads. ! TODO: Forces need to be appropriately scaled as well! do iw = 1, nbeads - eclas = eclas + Epot(iw) + Eclas = Eclas + Epot(iw) end do - eclas = eclas / nbeads + Eclas = Eclas / nbeads + + call numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads) + end subroutine force_h2o_schwenke ! TODO: Implement numerical forces generally for all potentials From dc12eb6fac1119c666dedd60d899279bef69c9a6 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Fri, 26 Jan 2024 19:04:12 +0000 Subject: [PATCH 14/15] Actually test --- tests/test.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test.sh b/tests/test.sh index e2679cac..b4f57a0b 100755 --- a/tests/test.sh +++ b/tests/test.sh @@ -66,7 +66,7 @@ function diff_files { if [[ $error_code -ne 0 ]];then # The reference file is different, but maybe it's just numerical noise? error_code=0 - diff -y -W 500 $test_file $ref_file | egrep -e '|' -e '<' -e '>' > $test_file.diff + diff -y -W 500 $test_file $ref_file | grep -e '|' -e '<' -e '>' > $test_file.diff ../numdiff.py $test_file.diff || error_code=$? @@ -127,7 +127,7 @@ if [[ $TESTS = "all" ]];then LZ_SS LZ_ST LZ_ENE \ PIMD ABINITIO ABINITIO-FAIL MTS \ LANGEVIN QT QT2 PIGLE PIGLE2 GLE-CANONICAL \ - HARMON MORSE DOUBLEWELL SPLINE MM MINI QMMM \ + H2O_SCHWENKE HARMON MORSE DOUBLEWELL SPLINE MM MINI QMMM \ ANALYZE_EXT CMDLINE WATER_FAIL ERMD) if [[ $MPI = "TRUE" ]];then From f3393388ccf81af444f644a9988c76422f821870 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Fri, 26 Jan 2024 19:19:47 +0000 Subject: [PATCH 15/15] Remove original file --- h2o_schwenke_orig.f | 493 -------------------------------------------- src/force_h2o.F90 | 9 +- 2 files changed, 6 insertions(+), 496 deletions(-) delete mode 100644 h2o_schwenke_orig.f diff --git a/h2o_schwenke_orig.f b/h2o_schwenke_orig.f deleted file mode 100644 index 57a4df53..00000000 --- a/h2o_schwenke_orig.f +++ /dev/null @@ -1,493 +0,0 @@ - program h2pes - implicit real*8 (a-h,o-z) - dimension rij(1,3),v(1) - rij(1,1) = 1.8 - rij(1,2) = 1.8 - rij(1,3) = 1.7 - call vibpot(rij, v, 1) - print *, "V = ", v(1) - contains - subroutine vibpot(rij,v,n) -c -c pes for h2o, -c Harry Partridge and David W. Schwenke, J. Chem. Phys., -c submitted Aug. 28, 1996. -c rij(i,1)& rij(i,2) are oh distances in au -c rij(i,3) is hoh angle in rad -c v(i) is pes in au -c n is number of geometries -c - dimension rij(n,3),v(n),c5z(245),cbasis(245),ccore(245), - $ crest(245),idx(245,3),fmat(15,3) - data (idx(i,1),i=1,245)/ - $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, - $ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, - $ 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, - $ 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, - $ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - $ 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, - $ 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, - $ 6, 6, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, - $ 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 5, 5, - $ 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, - $ 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, - $ 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 9, 9, - $ 9, 9, 9, 9, 9/ - data (idx(i,2),i=1,245)/ - $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, - $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, - $ 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, - $ 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, - $ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, - $ 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, - $ 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, - $ 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, - $ 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, - $ 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, - $ 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, - $ 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 1, 1, - $ 1, 1, 1, 1, 1/ - data (idx(i,3),i=1,245)/ - $ 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15, 1, 2, 3, 4, 5, - $ 6, 7, 8, 9,10,11,12,13,14, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11, - $12,13, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13, 1, 2, 3, 4, 5, - $ 6, 7, 8, 9,10,11,12, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12, 1, - $ 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, - $11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8, - $ 9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 1, 2, 3, 4, 5, 6, 7, 8, - $ 9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, - $ 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, - $ 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, - $ 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, - $ 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, - $ 3, 4, 5, 6, 7/ - data (c5z(i),i=1,245)/ - $ 4.2278462684916D+04, 4.5859382909906D-02, 9.4804986183058D+03, - $ 7.5485566680955D+02, 1.9865052511496D+03, 4.3768071560862D+02, - $ 1.4466054104131D+03, 1.3591924557890D+02,-1.4299027252645D+03, - $ 6.6966329416373D+02, 3.8065088734195D+03,-5.0582552618154D+02, - $-3.2067534385604D+03, 6.9673382568135D+02, 1.6789085874578D+03, - $-3.5387509130093D+03,-1.2902326455736D+04,-6.4271125232353D+03, - $-6.9346876863641D+03,-4.9765266152649D+02,-3.4380943579627D+03, - $ 3.9925274973255D+03,-1.2703668547457D+04,-1.5831591056092D+04, - $ 2.9431777405339D+04, 2.5071411925779D+04,-4.8518811956397D+04, - $-1.4430705306580D+04, 2.5844109323395D+04,-2.3371683301770D+03, - $ 1.2333872678202D+04, 6.6525207018832D+03,-2.0884209672231D+03, - $-6.3008463062877D+03, 4.2548148298119D+04, 2.1561445953347D+04, - $-1.5517277060400D+05, 2.9277086555691D+04, 2.6154026873478D+05, - $-1.3093666159230D+05,-1.6260425387088D+05, 1.2311652217133D+05, - $-5.1764697159603D+04, 2.5287599662992D+03, 3.0114701659513D+04, - $-2.0580084492150D+03, 3.3617940269402D+04, 1.3503379582016D+04, - $-1.0401149481887D+05,-6.3248258344140D+04, 2.4576697811922D+05, - $ 8.9685253338525D+04,-2.3910076031416D+05,-6.5265145723160D+04, - $ 8.9184290973880D+04,-8.0850272976101D+03,-3.1054961140464D+04, - $-1.3684354599285D+04, 9.3754012976495D+03,-7.4676475789329D+04, - $-1.8122270942076D+05, 2.6987309391410D+05, 4.0582251904706D+05, - $-4.7103517814752D+05,-3.6115503974010D+05, 3.2284775325099D+05, - $ 1.3264691929787D+04, 1.8025253924335D+05,-1.2235925565102D+04, - $-9.1363898120735D+03,-4.1294242946858D+04,-3.4995730900098D+04, - $ 3.1769893347165D+05, 2.8395605362570D+05,-1.0784536354219D+06, - $-5.9451106980882D+05, 1.5215430060937D+06, 4.5943167339298D+05, - $-7.9957883936866D+05,-9.2432840622294D+04, 5.5825423140341D+03, - $ 3.0673594098716D+03, 8.7439532014842D+04, 1.9113438435651D+05, - $-3.4306742659939D+05,-3.0711488132651D+05, 6.2118702580693D+05, - $-1.5805976377422D+04,-4.2038045404190D+05, 3.4847108834282D+05, - $-1.3486811106770D+04, 3.1256632170871D+04, 5.3344700235019D+03, - $ 2.6384242145376D+04, 1.2917121516510D+05,-1.3160848301195D+05, - $-4.5853998051192D+05, 3.5760105069089D+05, 6.4570143281747D+05, - $-3.6980075904167D+05,-3.2941029518332D+05,-3.5042507366553D+05, - $ 2.1513919629391D+03, 6.3403845616538D+04, 6.2152822008047D+04, - $-4.8805335375295D+05,-6.3261951398766D+05, 1.8433340786742D+06, - $ 1.4650263449690D+06,-2.9204939728308D+06,-1.1011338105757D+06, - $ 1.7270664922758D+06, 3.4925947462024D+05,-1.9526251371308D+04, - $-3.2271030511683D+04,-3.7601575719875D+05, 1.8295007005531D+05, - $ 1.5005699079799D+06,-1.2350076538617D+06,-1.8221938812193D+06, - $ 1.5438780841786D+06,-3.2729150692367D+03, 1.0546285883943D+04, - $-4.7118461673723D+04,-1.1458551385925D+05, 2.7704588008958D+05, - $ 7.4145816862032D+05,-6.6864945408289D+05,-1.6992324545166D+06, - $ 6.7487333473248D+05, 1.4361670430046D+06,-2.0837555267331D+05, - $ 4.7678355561019D+05,-1.5194821786066D+04,-1.1987249931134D+05, - $ 1.3007675671713D+05, 9.6641544907323D+05,-5.3379849922258D+05, - $-2.4303858824867D+06, 1.5261649025605D+06, 2.0186755858342D+06, - $-1.6429544469130D+06,-1.7921520714752D+04, 1.4125624734639D+04, - $-2.5345006031695D+04, 1.7853375909076D+05,-5.4318156343922D+04, - $-3.6889685715963D+05, 4.2449670705837D+05, 3.5020329799394D+05, - $ 9.3825886484788D+03,-8.0012127425648D+05, 9.8554789856472D+04, - $ 4.9210554266522D+05,-6.4038493953446D+05,-2.8398085766046D+06, - $ 2.1390360019254D+06, 6.3452935017176D+06,-2.3677386290925D+06, - $-3.9697874352050D+06,-1.9490691547041D+04, 4.4213579019433D+04, - $ 1.6113884156437D+05,-7.1247665213713D+05,-1.1808376404616D+06, - $ 3.0815171952564D+06, 1.3519809705593D+06,-3.4457898745450D+06, - $ 2.0705775494050D+05,-4.3778169926622D+05, 8.7041260169714D+03, - $ 1.8982512628535D+05,-2.9708215504578D+05,-8.8213012222074D+05, - $ 8.6031109049755D+05, 1.0968800857081D+06,-1.0114716732602D+06, - $ 1.9367263614108D+05, 2.8678295007137D+05,-9.4347729862989D+04, - $ 4.4154039394108D+04, 5.3686756196439D+05, 1.7254041770855D+05, - $-2.5310674462399D+06,-2.0381171865455D+06, 3.3780796258176D+06, - $ 7.8836220768478D+05,-1.5307728782887D+05,-3.7573362053757D+05, - $ 1.0124501604626D+06, 2.0929686545723D+06,-5.7305706586465D+06, - $-2.6200352535413D+06, 7.1543745536691D+06,-1.9733601879064D+04, - $ 8.5273008477607D+04, 6.1062454495045D+04,-2.2642508675984D+05, - $ 2.4581653864150D+05,-9.0376851105383D+05,-4.4367930945690D+05, - $ 1.5740351463593D+06, 2.4563041445249D+05,-3.4697646046367D+03, - $-2.1391370322552D+05, 4.2358948404842D+05, 5.6270081955003D+05, - $-8.5007851251980D+05,-6.1182429537130D+05, 5.6690751824341D+05, - $-3.5617502919487D+05,-8.1875263381402D+02,-2.4506258140060D+05, - $ 2.5830513731509D+05, 6.0646114465433D+05,-6.9676584616955D+05, - $ 5.1937406389690D+05, 1.7261913546007D+05,-1.7405787307472D+04, - $-3.8301842660567D+05, 5.4227693205154D+05, 2.5442083515211D+06, - $-1.1837755702370D+06,-1.9381959088092D+06,-4.0642141553575D+05, - $ 1.1840693827934D+04,-1.5334500255967D+05, 4.9098619510989D+05, - $ 6.1688992640977D+05, 2.2351144690009D+05,-1.8550462739570D+06, - $ 9.6815110649918D+03,-8.1526584681055D+04,-8.0810433155289D+04, - $ 3.4520506615177D+05, 2.5509863381419D+05,-1.3331224992157D+05, - $-4.3119301071653D+05,-5.9818343115856D+04, 1.7863692414573D+03, - $ 8.9440694919836D+04,-2.5558967650731D+05,-2.2130423988459D+04, - $ 4.4973674518316D+05,-2.2094939343618D+05/ - data (cbasis(i),i=1,245)/ - $ 6.9770019624764D-04,-2.4209870001642D+01, 1.8113927151562D+01, - $ 3.5107416275981D+01,-5.4600021126735D+00,-4.8731149608386D+01, - $ 3.6007189184766D+01, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $-7.7178474355102D+01,-3.8460795013977D+01,-4.6622480912340D+01, - $ 5.5684951167513D+01, 1.2274939911242D+02,-1.4325154752086D+02, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00,-6.0800589055949D+00, - $ 8.6171499453475D+01,-8.4066835441327D+01,-5.8228085624620D+01, - $ 2.0237393793875D+02, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 3.3525582670313D+02, 7.0056962392208D+01,-4.5312502936708D+01, - $-3.0441141194247D+02, 2.8111438108965D+02, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00,-1.2983583774779D+02, 3.9781671212935D+01, - $-6.6793945229609D+01,-1.9259805675433D+02, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00,-8.2855757669957D+02,-5.7003072730941D+01, - $-3.5604806670066D+01, 9.6277766002709D+01, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 8.8645622149112D+02,-7.6908409772041D+01, - $ 6.8111763314154D+01, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 2.5090493428062D+02,-2.3622141780572D+02, 5.8155647658455D+02, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 2.8919570295095D+03, - $-1.7871014635921D+02,-1.3515667622500D+02, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00,-3.6965613754734D+03, 2.1148158286617D+02, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00,-1.4795670139431D+03, - $ 3.6210798138768D+02, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $-5.3552886800881D+03, 3.1006384016202D+02, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 1.6241824368764D+03, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 4.3764909606382D+03, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 1.0940849243716D+03, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 3.0743267832931D+03, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00/ - data (ccore(i),i=1,245)/ - $ 2.4332191647159D-02,-2.9749090113656D+01, 1.8638980892831D+01, - $-6.1272361746520D+00, 2.1567487597605D+00,-1.5552044084945D+01, - $ 8.9752150543954D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $-3.5693557878741D+02,-3.0398393196894D+00,-6.5936553294576D+00, - $ 1.6056619388911D+01, 7.8061422868204D+01,-8.6270891686359D+01, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00,-3.1688002530217D+01, - $ 3.7586725583944D+01,-3.2725765966657D+01,-5.6458213299259D+00, - $ 2.1502613314595D+01, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 5.2789943583277D+02,-4.2461079404962D+00,-2.4937638543122D+01, - $-1.1963809321312D+02, 2.0240663228078D+02, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00,-6.2574211352272D+02,-6.9617539465382D+00, - $-5.9440243471241D+01, 1.4944220180218D+01, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00,-1.2851139918332D+03,-6.5043516710835D+00, - $ 4.0410829440249D+01,-6.7162452402027D+01, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 1.0031942127832D+03, 7.6137226541944D+01, - $-2.7279242226902D+01, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $-3.3059000871075D+01, 2.4384498749480D+01,-1.4597931874215D+02, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 1.6559579606045D+03, - $ 1.5038996611400D+02,-7.3865347730818D+01, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00,-1.9738401290808D+03,-1.4149993809415D+02, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00,-1.2756627454888D+02, - $ 4.1487702227579D+01, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $-1.7406770966429D+03,-9.3812204399266D+01, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00,-1.1890301282216D+03, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 2.3723447727360D+03, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00,-1.0279968223292D+03, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 5.7153838472603D+02, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00/ - data (crest(i),i=1,245)/ - $ 0.0000000000000D+00,-4.7430930170000D+00,-1.4422132560000D+01, - $-1.8061146510000D+01, 7.5186735000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $-2.7962099800000D+02, 1.7616414260000D+01,-9.9741392630000D+01, - $ 7.1402447000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00,-7.8571336480000D+01, - $ 5.2434353250000D+01, 7.7696745000000D+01, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 1.7799123760000D+02, 1.4564532380000D+02, 2.2347226000000D+02, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00,-4.3823284100000D+02,-7.2846553000000D+02, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00,-2.6752313750000D+02, 3.6170310000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, - $ 0.0000000000000D+00, 0.0000000000000D+00/ - data reoh,thetae,b1,roh,alphaoh,deoh,phh1,phh2/0.958649d0, - $ 104.3475d0,2.0d0,0.9519607159623009d0,2.587949757553683d0, - $ 42290.92019288289d0,16.94879431193463d0,12.66426998162947d0/ - data f5z,fbasis,fcore,frest/0.99967788500000d0, - $ 0.15860145369897d0,-1.6351695982132d0,1d0/ - save - data ifirst/0/ - if(ifirst.eq.0)then - ifirst=1 - write(6,1) - 1 format(/1x,'pes for h2o', - $ /1x,'by Harry Partridge and David W. Schwenke', - $ /1x,'submitted to J. Chem. Phys. Aug. 28, 1996') - write(6,56) - 56 format(/1x,'parameters before adjustment') - write(6,55)phh1,phh2,deoh,alphaoh,roh - 55 format(/1x,'two body potential parameters:', - $ /1x,'hh: phh1 = ',f10.1,' phh2 = ',f5.2, - $ /1x,'oh: deoh = ',f10.1,' alpha = ',f7.4, - $ ' re = ',f7.4) - write(6,4)reoh,thetae,b1 - 4 format(/1x,'three body parameters:', - $ /1x,'reoh = ',f10.4,' thetae = ',f10.4, - $ /1x,'betaoh = ',f10.4, - $ /1x,' i j k',7x,'c5z',9x,'cbasis',10x,'ccore', - $ 10x,'crest') - do 2 i=1,245 - write(6,5)(idx(i,j)-1,j=1,3),c5z(i),cbasis(i),ccore(i),crest(i) - 5 format(1x,3i5,1p4e15.7) - c5z(i)=f5z*c5z(i)+fbasis*cbasis(i)+fcore*ccore(i) - $ +frest*crest(i) - 2 continue - write(6,57)f5z,fbasis,fcore,frest - 57 format(/1x,'adjusting parameters using scale factors ', - $ /1x,'f5z = ',f11.8, - $ /1x,'fbasis = ',f11.8, - $ /1x,'fcore = ',f11.8, - $ /1x,'frest = ',f11.8) - phh1=phh1*f5z - deoh=deoh*f5z - write(6,55)phh1,phh2,deoh,alphaoh,roh - write(6,58)reoh,thetae,b1,((idx(i,j)-1,j=1,3),c5z(i),i=1,245) - 58 format(/1x,'three body parameters:', - $ /1x,'reoh = ',f10.4,' thetae = ',f10.4, - $ /1x,'betaoh = ',f10.4, - $ /1x,' i j k cijk', - $ /(1x,3i5,1pe15.7)) -c -c convert parameters from 1/cm, angstrom to a.u. -c - reoh=reoh/0.529177249d0 - b1=b1*0.529177249d0*0.529177249d0 - do 3 i=1,245 - c5z(i)=c5z(i)*4.556335d-6 - 3 continue - rad=acos(-1d0)/1.8d2 - ce=cos(thetae*rad) - phh1=phh1*exp(phh2) - phh1=phh1*4.556335d-6 - phh2=phh2*0.529177249d0 - deoh=deoh*4.556335d-6 - roh=roh/0.529177249d0 - alphaoh=alphaoh*0.529177249d0 - c5z(1)=c5z(1)*2d0 - end if - do 6 i=1,n - x1=(rij(i,1)-reoh)/reoh - x2=(rij(i,2)-reoh)/reoh - x3=cos(rij(i,3))-ce - rhh=sqrt(rij(i,1)**2+rij(i,2)**2 - $ -2d0*rij(i,1)*rij(i,2)*cos(rij(i,3))) - vhh=phh1*exp(-phh2*rhh) - ex=exp(-alphaoh*(rij(i,1)-roh)) - voh1=deoh*ex*(ex-2d0) - ex=exp(-alphaoh*(rij(i,2)-roh)) - voh2=deoh*ex*(ex-2d0) - fmat(1,1)=1d0 - fmat(1,2)=1d0 - fmat(1,3)=1d0 - do 10 j=2,15 - fmat(j,1)=fmat(j-1,1)*x1 - fmat(j,2)=fmat(j-1,2)*x2 - fmat(j,3)=fmat(j-1,3)*x3 - 10 continue - v(i)=0d0 - do 12 j=2,245 - term=c5z(j)*(fmat(idx(j,1),1)*fmat(idx(j,2),2) - $ +fmat(idx(j,2),1)*fmat(idx(j,1),2)) - $ *fmat(idx(j,3),3) - v(i)=v(i)+term - 12 continue - v(i)=v(i)*exp(-b1*((rij(i,1)-reoh)**2+(rij(i,2)-reoh)**2)) - $ +c5z(1) - $ +voh1+voh2+vhh - 6 continue - return - end - end diff --git a/src/force_h2o.F90 b/src/force_h2o.F90 index 1d7c3667..a9305b9c 100644 --- a/src/force_h2o.F90 +++ b/src/force_h2o.F90 @@ -1,6 +1,6 @@ ! Interface to analytical H2O (single water molecule) potentials. ! This is invoked via pot='_h2o_', and the potential is selected via -! h2opot='schwenke' +! h2opot='schwenke' in namelist General in ABIN input file. module mod_force_h2o use mod_const, only: DP use mod_files, only: stdout, stderr @@ -27,10 +27,11 @@ subroutine initialize_h2o_pot(natom, atom_names) write (stderr, *) 'Water atoms must be ordered as "O H H"' call fatal_error(__FILE__, __LINE__, "This is not a water molecule!") end if + ! TODO: Some PES initialization could be performed here end subroutine initialize_h2o_pot - ! Compute internal coordinates from cartesians - ! OH distances in bohrs, HOH angle in radians + ! Compute internal coordinates from cartesians, + ! OH distances in Bohrs, HOH angle in radians. subroutine get_internal_coords(x, y, z, iw, rOH1, rOH2, aHOH_rad) use mod_const, only: PI use mod_utils, only: get_distance, get_angle @@ -81,7 +82,9 @@ subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, Eclas, natom, nbeads) end do ! Calculated the potential for all PI beads at once + ! The potential is implemented in h2o_schwenke.f call h2o_pot_schwenke(rij, Epot, nbeads) + ! For Path Integrals, the final energy of the PI necklace ! is an average over all beads. ! TODO: Forces need to be appropriately scaled as well!