diff --git a/.fprettify.rc b/.fprettify.rc index 7ad662da..8cc6fd9e 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, h2o_schwenke.f] +exclude=[random.F90,fftw3.F90,force_cp2k.F90, h2o_schwenke.f, h2o_cvrqd.f] diff --git a/src/Makefile b/src/Makefile index ff25d084..c1d6422b 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 h2o_schwenke.o force_h2o.o force_cp2k.o sh_integ.o surfacehop.o landau_zener.o\ + remd.o force_bound.o water.o h2o_schwenke.o h2o_cvrqd.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 diff --git a/src/force_h2o.F90 b/src/force_h2o.F90 index a9305b9c..bf5d945b 100644 --- a/src/force_h2o.F90 +++ b/src/force_h2o.F90 @@ -57,6 +57,9 @@ subroutine force_h2o(x, y, z, fx, fy, fz, eclas, natom, nbeads) ! 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 if (h2opot == 'cvrqd') then + call force_h2o_cvrqd(x, y, z, fx, fy, fz, eclas, natom, nbeads) + call fatal_error(__FILE__, __LINE__, 'Numerical forces not yet implemented!') else call fatal_error(__FILE__, __LINE__, 'Potential '//trim(h2opot)//' not implemented') end if @@ -97,6 +100,39 @@ subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, Eclas, natom, nbeads) end subroutine force_h2o_schwenke + subroutine force_h2o_cvrqd(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) :: E + integer :: iw + real(DP) :: mH, mO + + ! TODO: Pass in the actual masses + ! What should be these set to? Should we set pure isotopes? + mH = 1.008D0 + mO = 15.999D0 + + Eclas = 0.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) + + call h2o_pot_cvrqd(E, rOH1, rOH2, aHOH_rad, mO, mH) + + Eclas = Eclas + E + end do + Eclas = Eclas / nbeads + + ! TODO: Given the small difference between the Schwenke potential, + ! we might not need to implement numerical forces here. + ! call numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads) + + end subroutine force_h2o_cvrqd + ! 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) diff --git a/src/fortran_interfaces.F90 b/src/fortran_interfaces.F90 index 696eb911..5648a7e2 100644 --- a/src/fortran_interfaces.F90 +++ b/src/fortran_interfaces.F90 @@ -63,14 +63,21 @@ function usleep(useconds) bind(c, name='usleep') integer(kind=C_INT) :: usleep end function usleep - ! Returns a potential energy of a water molecule + ! Computes potential energy of a water molecule ! using Schwenke potential, see h2o_schwenke.f subroutine h2o_pot_schwenke(rij, v, n) import :: DP integer, intent(in) :: n - real(DP) :: rij(n, 3), v(n) + real(DP), intent(in) :: rij(n, 3) + real(DP), intent(out) :: v(n) end subroutine h2o_pot_schwenke + subroutine h2o_pot_cvrqd(V, rOH1, rOH2, aHOH, mH, mO) + import :: DP + real(DP), intent(out) :: V + real(DP), intent(in) :: rOH1, rOH2, aHOH, mH, mO + end subroutine h2o_pot_cvrqd + end interface end module mod_interfaces diff --git a/src/h2o_cvrqd.f b/src/h2o_cvrqd.f new file mode 100644 index 00000000..b1e24052 --- /dev/null +++ b/src/h2o_cvrqd.f @@ -0,0 +1,1544 @@ +! Single Water PES published in: +! CVRQD ab initio ground-state adiabatic potential energy surfaces for the water molecule +! https://doi.org/10.1063/1.2378766 + +! inputs are in bohr for the OH distances Q1,Q2, +! and in radians for the H-O-H angle . +! The parameters xma and xmb are the mass of nucleus of O and H, +! respectively. They are read from the dvr input. + SUBROUTINE h2o_pot_cvrqd(V,Q1,Q2,theta,xma,xmb) + implicit double precision(a-h,o-y),logical(z) + + CALL PESleq6(VQ,Q1,Q2,THETA) + CALL BREITB3lin(Vbr,Q1,Q2,THETA) + CALL pesd2x(Vd,Q1,Q2,THETA) + CALL bodc16(V16,Q1,Q2,THETA) + CALL bodc18(V18,Q1,Q2,THETA) + xmaso16=15.990526d0 + xmaso17=16.9947425d0 + xmaso18=17.9947714d0 + xmasd=2.013553214d0 + xmash=1.00727647d0 + v1=(xmaso18*xmaso16*(v16-v18)/(xmaso18-xmaso16))/xma+ + *(2*xmash*(xmaso18*v18-xmaso16*v16)/(2*(xmaso18-xmaso16)))/xmb + + CALL rel(Vr,Q1,Q2,THETA) + att1=Q1*0.5291772d0 + att2=Q2*0.5291772d0 + CALL wifin(att1,att2,THETA,vp) + CALL cvps(Vcvps,Q1,Q2,THETA) + +! This is where the ab initio potential gets defined from the bits and pieces + + v=vp/219474.624d0+Vcvps+vr+vbr+vd+VQ+V1 + +! vp is the extrapolated ICMRCI basic potential +! vr is the relativistic (MVD1) incremental potential +! vbr is the Breit increment +! vd is the MVD2 increment +! v1 is the BODC increment +! vq is the QED increment + + RETURN + END SUBROUTINE h2o_pot_cvrqd + + SUBROUTINE BODC16(V,R1,R2,THETA) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION FT(69), cv(69) + integer :: i, i5 + integer :: nv, ZERO + + DATA C0/0.0D0/,SCALE/1.0D-6/ + data nv/69/ + DATA ZERO/0.0D0/ + DATA RZ/0.9576257 /,RHO/75.48992362 / + DATA TOANG/0.5291772/ + DATA cv/2785.4109260806,-90.0463360612,48.5226385100, + *201.6994168988,-19.4625405248,217.8639531832, + *-48.9272565544,-76.0803927294,25.9675092231,-65.2206704082, + *60.7284555026,-386.5173198637,-2.0352745279,3.5405842322, + *-32.1830475213,126.8885218953,60.2606780157,134.6300248947, + *-14.0991224631,560.8792812623,27.8464660146,-67.1464124186, + *-181.5695164733,43.7287769150,-150.7468776877,-103.2652779984, + *-159.6432245839,-24.8956181482,185.4825335893,-231.4775497546, + *-51.9208759775,3.7066140140,-212.4422129941,207.1884490703, + *383.1659239100,50.8641728660,35.1754939408,127.2280510484, + *-154.4544312699,55.5787967758,282.6216945851,116.5606405651, + *5.4433254144,-107.1461094167,-173.8895717556,-26.5859674990, + *-560.8756697840,-237.7109212157,143.9462552048,-592.3478209334, + *0.0,-198.6468835005,-19.9674473372,-14.1731270087,193.4510720304, + *4.6347021028,32.9502486772,-221.1685318724,26.4090449111, + *-268.432837385,-147.1422366151,133.5465868568,363.9554096142, + *673.9006484856,214.9454642643,40.7735822438,65.2246188257, + *173.0708970426,1.9795259929/ + + THETAeq = (180.d0 - RHO)*3.141592654/180.d0 + Req= RZ/toang + S1 = R1 - Req + S2 = DCOS(THETA) - DCOS(THETAeq) + S3 = R2 - Req + + Y1 = (S1 + S3)/2.0d0 + Y2 = S2 + Y3 = (S1 - S3)/2.0d0 + + DO I5=1,NV + FT(I5)=0.0D0 + END DO + + FT(1)=1.0d0 + + FT(2)=Y1 + FT(3)=Y2 + + FT(4)=Y1**2 + FT(5)=Y2**2 + FT(6)=Y3**2 + FT(7)=Y1*Y2 + + FT(8)=Y1**3 + FT(9)=Y2**3 + FT(10)=Y1**2*Y2 + FT(11)=Y2**2*Y1 + FT(12)=Y3**2*Y1 + FT(13)=Y3**2*Y2 + + FT(14)=Y1**4 + FT(15)=Y2**4 + FT(16)=Y3**4 + FT(17)=Y1**3*Y2 + FT(18)=Y2**3*Y1 + FT(19)=Y1**2*Y2**2 + FT(20)=Y1**2*Y3**2 + FT(21)=Y2**2*Y3**2 + FT(22)=Y3**2*Y1*Y2 + + FT(23)=Y1**5 + FT(24)=Y2**5 + FT(25)=Y1**4*Y2 + FT(26)=Y2**4*Y1 + FT(27)=Y3**4*Y1 + FT(28)=Y3**4*Y2 + FT(29)=Y1**3*Y2**2 + FT(30)=Y1**3*Y3**2 + FT(31)=Y2**3*Y1**2 + FT(32)=Y2**3*Y3**2 + FT(33)=Y1**2*Y2*Y3**2 + FT(34)=Y1*Y2**2*Y3**2 + + FT(35)=Y1**6 + FT(36)=Y2**6 + FT(37)=Y3**6 + FT(38)=Y1**5*Y2 + FT(39)=Y2**5*Y1 + FT(40)=Y1**4*Y2**2 + FT(41)=Y2**4*Y1**2 + FT(42)=Y2**4*Y3**2 + FT(43)=Y3**4*Y2**2 + FT(44)=Y1**4*Y3**2 + FT(45)=Y3**4*Y1**2 + FT(46)=Y3**4*Y1*Y2 + FT(47)=Y1**3*Y2**3 + FT(48)=Y2**3*Y1**2*Y3**2 + FT(49)=Y1**3*Y3**2*Y2 + FT(50)=Y2**3*Y3**2*Y1 + FT(51)=Y1**2*Y2**2*Y1**2 + + + FT(52)=Y1**7 + FT(53)=Y2**7 + FT(54)=Y1**6*Y2 + FT(55)=Y2**6*Y1 + FT(56)=Y3**6*Y1 + FT(57)=Y3**6*Y2 + FT(58)=Y1**5*Y2**2 + FT(59)=Y1**5*Y3**2 + FT(60)=Y2**5*Y1**2 + FT(61)=Y2**5*Y3**2 + FT(62)=Y1**4*Y2*Y3**2 + FT(63)=Y1**4*Y2**3 + FT(64)=Y2**4*Y1*Y3**2 + FT(65)=Y2**4*Y1**3 + FT(66)=Y3**4*Y1*Y2**2 + FT(67)=Y3**4*Y2*Y1**2 + FT(68)=Y3**4*Y1**3 + FT(69)=Y3**4*Y2**3 + + V=ZERO + DO I=1,NV + V=V+CV(I)*FT(I) + END DO + V=C0+SCALE*V + RETURN + END SUBROUTINE bodc16 + + SUBROUTINE BODC18(V,R1,R2,THETA) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION FT(117), cv(117) + integer :: nv, i, i5 + + DATA C0/0.0D0/,SCALE/1.0D-6/ + data nv/117/ + DATA ZERO/0.0D0/ + DATA RZ/0.9576257 /,RHO/75.48992362 / + DATA TOANG/0.5291772/ + DATA cv/2520.8349350538,-91.1209324309,49.9198880161, + *201.5986877584,-27.0761946057,213.0632778654,-42.3461775982, + *-89.5081758955,3.1590606555,-82.9610483797,91.2462747271, + *-345.6845385214,2.9125633339,-9.2743637228,19.1249589270, + *126.8439173732,102.4748093797,29.1206937761,-0.4903235370, + *548.5647508851,57.4425043137,-65.3115364766,-32.7831814238, + *111.9709105628,-278.4048997815,-231.4930911427,-323.4608284981, + *-98.0602324965,181.1676886771,-282.2900662875,167.5274873935, + *47.8266275957,-213.7402930366,20.8555467237,169.6474418937, + *-62.1263131255,71.9040085206,72.1511749317,123.9432784777, + *173.6777915640,58.5107154936,-27.8257872910,99.6002434065, + *-135.2045098416,145.2475743775,-23.3633794048,-599.5375919902, + *-540.4989859645,275.5466454040,-502.8009780288,167.2307803194, + *-109.3097362364,-77.1423552840,492.0803282409,313.6574722992, + *-22.4885664677,141.5175472800,-1134.9499411146,88.6368503668, + *-587.9451428738,-247.5329792460,-4.9333140627,322.2004639448, + *916.1768044740,479.0204491692,93.9437933859,107.5055433254, + *0.0,63.1640462956,191.6712842623,1.6032078252,78.8696064123, + *1.0174252672,-377.6585243760,-192.1212462644,653.6712620983, + *0.0,261.5131272695,146.4973633837,0.0,-137.1867814770, + *-98.7015480745,0.0,0.0,0.0,0.0,-97.1438825386,0.0,0.0,0.0,0.0, + *0.0,22.7027546667,0.0,0.0,0.0,-341.9933632553,0.0,0.0, + *-111.8520186501,0.0,0.0,0.0,0.0,-48.2716070600,0.0,0.0,0.0,0.0, + *-5.0447226741,0.0,0.0,0.0,0.0,0.0,0.0,144.1528855579/ + THETAeq = (180.d0 - RHO)*3.141592654/180.d0 + Req= RZ/toang + S1 = R1 - Req + S2 = DCOS(THETA) - DCOS(THETAeq) + S3 = R2 - Req + + Y1 = (S1 + S3)/2.0d0 + Y2 = S2 + Y3 = (S1 - S3)/2.0d0 + +c Potential in Y1,Y2,Y3. Full set of +c derivatives for each order (1-4) and upper terms for others. + + DO I5=1,NV + FT(I5)=0.0D0 + END DO + + FT(1)=1.0d0 + + FT(2)=Y1 + FT(3)=Y2 + + FT(4)=Y1**2 + FT(5)=Y2**2 + FT(6)=Y3**2 + FT(7)=Y1*Y2 + + FT(8)=Y1**3 + FT(9)=Y2**3 + FT(10)=Y1**2*Y2 + FT(11)=Y2**2*Y1 + FT(12)=Y3**2*Y1 + FT(13)=Y3**2*Y2 + + FT(14)=Y1**4 + FT(15)=Y2**4 + FT(16)=Y3**4 + FT(17)=Y1**3*Y2 + FT(18)=Y2**3*Y1 + FT(19)=Y1**2*Y2**2 + FT(20)=Y1**2*Y3**2 + FT(21)=Y2**2*Y3**2 + FT(22)=Y3**2*Y1*Y2 + + FT(23)=Y1**5 + FT(24)=Y2**5 + FT(25)=Y1**4*Y2 + FT(26)=Y2**4*Y1 + FT(27)=Y3**4*Y1 + FT(28)=Y3**4*Y2 + FT(29)=Y1**3*Y2**2 + FT(30)=Y1**3*Y3**2 + FT(31)=Y2**3*Y1**2 + FT(32)=Y2**3*Y3**2 + FT(33)=Y1**2*Y2*Y3**2 + FT(34)=Y1*Y2**2*Y3**2 + + FT(35)=Y1**6 + FT(36)=Y2**6 + FT(37)=Y3**6 + FT(38)=Y1**5*Y2 + FT(39)=Y2**5*Y1 + FT(40)=Y1**4*Y2**2 + FT(41)=Y2**4*Y1**2 + FT(42)=Y2**4*Y3**2 + FT(43)=Y3**4*Y2**2 + FT(44)=Y1**4*Y3**2 + FT(45)=Y3**4*Y1**2 + FT(46)=Y3**4*Y1*Y2 + FT(47)=Y1**3*Y2**3 + FT(48)=Y2**3*Y1**2*Y3**2 + FT(49)=Y1**3*Y3**2*Y2 + FT(50)=Y2**3*Y3**2*Y1 + FT(51)=Y1**2*Y2**2*Y1**2 + + + FT(52)=Y1**7 + FT(53)=Y2**7 + FT(54)=Y1**6*Y2 + FT(55)=Y2**6*Y1 + FT(56)=Y3**6*Y1 + FT(57)=Y3**6*Y2 + FT(58)=Y1**5*Y2**2 + FT(59)=Y1**5*Y3**2 + FT(60)=Y2**5*Y1**2 + FT(61)=Y2**5*Y3**2 + FT(62)=Y1**4*Y2*Y3**2 + FT(63)=Y1**4*Y2**3 + FT(64)=Y2**4*Y1*Y3**2 + FT(65)=Y2**4*Y1**3 + FT(66)=Y3**4*Y1*Y2**2 + FT(67)=Y3**4*Y2*Y1**2 + FT(68)=Y3**4*Y1**3 + FT(69)=Y3**4*Y2**3 + FT(70)=Y1**3*Y3**2*Y2**2 + + FT(71)=Y1**8 + FT(72)=Y2**8 + FT(73)=Y3**8 + FT(74)=Y1**7*Y2 + FT(75)=Y2**7*Y1 + FT(76)=Y1**6*Y2**2 + FT(77)=Y1**6*Y3**2 + FT(78)=Y2**6*Y1**2 + FT(79)=Y2**6*Y3**2 + FT(80)=Y3**6*Y1**2 + FT(81)=Y3**6*Y2**2 + FT(82)=Y3**6*Y2*Y1 + FT(83)=Y1**5*Y3**2*Y2 + FT(84)=Y1**5*Y2**3 + FT(85)=Y2**5*Y1**3 + FT(86)=Y2**5*Y3**2*Y1 + FT(87)=Y1**4*Y2**4 + FT(88)=Y1**4*Y3**4 + FT(89)=Y1**4*Y2**2*Y3**2 + FT(90)=Y2**4*Y3**4 + FT(91)=Y2**4*Y3**2*Y1**2 + FT(92)=Y3**4*Y1**2*Y2**2 + FT(93)=Y3**4*Y1*Y2**3 + FT(94)=Y3**4*Y1**3*Y2 + FT(95)=Y1**3*Y2**3*Y3**2 + + FT(96)=Y2**9 + FT(97)=Y2**9*Y1*Y3**4 + FT(98)=Y2**9*Y1**5 + + FT(99)=Y2**10 + FT(100)=Y2**10*Y1**2*Y3**2 + FT(101)=Y2**10*Y3**4 + FT(102)=Y2**10*Y1**4 + FT(103)=Y2**11 + FT(104)=Y2**11*Y1*Y3**2 + FT(105)=Y2**11*Y1**3 + + FT(106)=Y2**12 + FT(107)=Y2**12*Y1**2 + FT(108)=Y2**12*Y3**2 + + FT(109)=Y2**13 + FT(110)=Y2**13*Y1 + + FT(111)=Y2**14 + + FT(112)=Y2**8*Y1**6 + FT(113)=Y2**8*Y3**6 + FT(114)=Y2**8*Y1**4*Y3**2 + FT(115)=Y2**8*Y3**4*Y1**2 + + FT(116)=Y2**7*Y1**7 + FT(117)=Y2**7*Y1**5*Y3**2 + + V=ZERO + DO I=1,NV + V=V+CV(I)*FT(I) + END DO +C SCALE AND SHIFT THE ZERO + V=C0+SCALE*V + RETURN + END + + SUBROUTINE rel(V,R1,R2,THETA) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION FT(82), cv(82) + integer :: nv, ZERO, i, i5, idx + + data nv/82/ + DATA ZERO/0.0D0/ + DATA RZ/0.9576257 /,RHO/75.48992362 / + DATA TOANG/0.5291772/ + DATA cv/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-0.0519929240, + *-0.0000574005,-0.0001860844,-0.0006152872,0.0000879624, + *-0.0004932809,-0.0000565213,0.0006886853,-0.0001027470, + *0.0003081795,0.0001078370,0.0019925277,0.0001286772, + *-0.0004511478,0.0000072157,-0.0004439080,-0.0001416903, + *0.0000984644,-0.0003092456,-0.0029118771,-0.0000562847, + *0.0000377062,0.0002610717,0.0000367337,0.0,-0.0002565517, + *0.0008262448,0.0,-0.0001036317,0.0025090209,0.0002105614, + *-0.0000204982,0.0,-0.0000414402,-0.0000676532,0.0000918240, + *-0.0000400594,0.0,-0.0006468368,0.0001397619,0.0005356017, + *0.0001585601,-0.0001899817,-0.0015914516,-0.0002822918, + *0.0,0.0004567782,0.0,0.0001064879,0.0,0.0,0.0,-0.0001250236, + *0.0000077559,0.0007064063,0.0,0.0,0.0,0.0006102666,-0.0004987995, + *0.0,-0.0001943536,-0.0002855510,0.0,-0.0002176976,0.0002410759, + *-0.0001644075,-0.0001182853,0.0,0.0,-0.0000272857,0.0/ + + THETAeq = (180.d0 - RHO)*3.141592654/180.d0 + Req= RZ/toang + S1 = R1 - Req + S2 = DCOS(THETA) - DCOS(THETAeq) + S3 = R2 - Req + + Y1 = (S1 + S3)/2.0d0 + Y2 = S2 + Y3 = (S1 - S3)/2.0d0 + + DO 88 I5=1,NV + FT(I5)=0.0D0 + 88 continue + FT(11)=1.0d0 + FT(12)=Y1 + FT(13)=Y2 + FT(14)=Y1**2 + FT(15)=Y2**2 + FT(16)=Y3**2 + FT(17)=Y1*Y2 + FT(18)=Y1**3 + FT(19)=Y2**3 + FT(20)=Y1**2*Y2 + FT(21)=Y2**2*Y1 + FT(22)=Y3**2*Y1 + FT(23)=Y3**2*Y2 + FT(24)=Y1**4 + FT(25)=Y2**4 + FT(26)=Y3**4 + FT(27)=Y1**3*Y2 + FT(28)=Y2**3*Y1 + FT(29)=Y1**2*Y2**2 + FT(30)=Y1**2*Y3**2 + FT(31)=Y2**2*Y3**2 + FT(32)=Y3**2*Y1*Y2 + FT(33)=Y1**5 + FT(34)=Y2**5 + FT(35)=Y1**4*Y2 + FT(36)=Y2**4*Y1 + FT(37)=Y3**4*Y1 + FT(38)=Y3**4*Y2 + FT(39)=Y1**3*Y2**2 + FT(40)=Y1**3*Y3**2 + FT(41)=Y2**3*Y1**2 + FT(42)=Y2**3*Y3**2 + FT(43)=Y1**2*Y2*Y3**2 + FT(44)=Y1*Y2**2*Y3**2 + FT(45)=Y1**6 + FT(46)=Y2**6 + FT(47)=Y3**6 + FT(48)=Y1**5*Y2 + FT(49)=Y2**5*Y1 + FT(50)=Y1**4*Y2**2 + FT(51)=Y2**4*Y1**2 + FT(52)=Y2**4*Y3**2 + FT(53)=Y3**4*Y2**2 + FT(54)=Y1**4*Y3**2 + FT(55)=Y3**4*Y1**2 + FT(56)=Y3**4*Y1*Y2 + FT(57)=Y1**3*Y2**3 + FT(58)=Y2**3*Y1**3 + FT(59)=Y1**3*Y3**2*Y2 + FT(60)=Y2**3*Y3**2*Y1 + FT(61)=Y1**2*Y2**2*Y1**2 + FT(62)=Y1**7 + FT(63)=Y2**7 + FT(64)=Y1**6*Y2 + FT(65)=Y2**6*Y1 + FT(66)=Y3**6*Y1 + FT(67)=Y3**6*Y2 + FT(68)=Y1**5*Y2**2 + FT(69)=Y1**5*Y3**2 + FT(70)=Y2**5*Y1**2 + FT(71)=Y2**5*Y3**2 + FT(72)=Y1**4*Y2*Y3**2 + FT(73)=Y1**4*Y2**3 + FT(74)=Y2**4*Y1*Y3**2 + FT(75)=Y2**4*Y1**3 + FT(76)=Y3**4*Y1*Y2**2 + FT(77)=Y3**4*Y2*Y1**2 + FT(78)=Y3**4*Y1**3 + FT(79)=Y3**4*Y2**3 + FT(80)=Y1**3*Y3**2*Y2**2 + FT(81)=Y1**8 + FT(82)=Y2**8 + + V=ZERO + DO 40 I=12,NV + V=V+CV(I)*FT(I) + 40 continue + RETURN + END + + SUBROUTINE BREITB3lin(VR,X,Y,Z) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + DATA SCALE/1.0D-6/ +c This subroutine contains corrections to the water NBO PES due to the BREIT +c term. +c [see for istance HM Quiney et al Chem. Phys. Lett. 290 (1998) 473 , +c Bethe and Salpheter, "Quantum mechanics of one and two-electron atoms"] +c The corrections have been computed on a grid based on the 325 point grid +c from P&S (see J. chem. phys., 106 (1997) 4618). +c Moreover, a few extra points have been added , as well as a cut in the radial +c coordinates (lines 1-2), in order to account for high bending modes. +c The final grid used contains 293 points. +c Then the points have been fitted with a polynomial in X,Y and Z, using a +c Mathematica script. +c The basis set used for the electronic calculations is the set called 'B' +c provided by H.M. Quiney (see Chem. Phys. Lett. ). +c +c Those corrections have been computed by Paolo +c email: paolo@theory.phys.ucl.ac.uk . + +c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +c The inputs are in u.a., and X,Y are the distances of the H atoms from +c the oxygen, and Z is the angle HOH in radiants. The final result is in +c Hartree. + +c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +c Potential symmetrical in X,Y,Z. Full set of +c derivatives for each order (1-4) and upper terms for others + + Seq=1.8240445d0 + e=1.809645d0 + h=0.5d0 + + if (x.gt.3.5d0) x=3.5d0 + if (y.gt.3.5d0) y=3.5d0 + + xm=x-y + xp=x+y + xp2=-e+h*xp + zm=-Seq+z + + v1=7577.371974812879+54.16425693884002*xm**2+ + y 6.774810680339735*xm**4-0.7125792752605449*xm**6- + y 212.3502737780579*xp2- + y 119.8400228527018*xm**2*xp2-29.18501086262345*xm**4*xp2+ + y 221.1271759652365*xp2**2+72.97844564936518*xm**2*xp2**2+ + y 89.0829000401668*xm**4*xp2**2-180.6736318109525*xp2**3+ + y 321.1679099895901*xm**2*xp2**3-86.9772033725647*xm**4*xp2**3+ + y 81.467133012329*xp2**4-890.357807854604*xm**2*xp2**4+ + y 57.52112099193129*xp2**5+676.0925261740228*xm**2*xp2**5- + y 58.228375665496*xp2**6+1.323061103513093*zm- + y 0.5270121336782832*xm**2*zm+0.01985434592556156*xp2* + y zm+7.715533877300036*xm**2*xp2*zm-6.229144779356836*xm**4*xp2* + y zm-0.1710351457311522*xp2**2*zm+9.92240371274571*xm**4*xp2**2* + y zm+10.47917821735698*xp2**3*zm-27.67291046310705*xm**2*xp2**4* + y zm+15.42094531691969*xp2**6*zm+7.749223096520598*zm**2- + y 0.2147443261198608*xm**2*zm**2-1.164069372530965*xm**4*zm**2+ + y 10.24148015693046*xp2*zm**2-2.364801830233726*xm**2*xp2*zm**2 + v2=3.405435470980123*xm**4*xp2*zm**2+11.53659470986242*xp2**2* + y zm**2+40.71096970108562*xp2**3*zm**2-65.49114275587444*xp2**4* + y zm**2+0.5246601257334035*zm**3+1.025008298074623*xm**2*zm**3+ + y 13.57824254391274*xp2*zm**3-7.469419914001589*xp2**2* + y zm**3-33.70757112970705*xp2**3*zm**3+30.20514216972833*xp2**4* + y zm**3-10.53913543447923*zm**4-0.6159136295163627*xm**2*zm**4- + y 19.56431274355461*xp2*zm**4-20.81965238209867*xp2**2*zm**4+ + y 5.998958874987758*xp2**3*zm**4-9.44711265431818*zm**5- + y 22.55622148750276*xp2*zm**5+16.30440168684215*xp2**2* + y zm**5+19.20675957512514*zm**6+19.78080962673524*xp2* + y zm**6+8.08923849773384*zm**7-10.68490632273025*zm**8 + + VR=V1+V2 + +C SCALE AND SHIFT THE ZERO + VR=VR*scale + + RETURN + END + + SUBROUTINE PESd2x(VR,x,y,z) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + +c This subroutine contains corrections to the water NBO PES due to the Darwin +c 2 electrons term. Those corrections have been computed by Gyorgy +c tarczay@para.chem.elte.hu +c Those corrections are computed on the P&S grid of 325 points. +c (see J. chem. phys., 106 (1997) 4618), to which a few x points have been +c added in irder to account to high bending modes. The final grid containd 341 +c points. + +c The input are in u.a., and X,Y are the distances of the H atoms from +c the oxygen, and Z is the angle HOH in radiants. The final result is in +c Hartree. + + Seq=1.8240445d0 + pi = dacos(-1.0d0) + e=1.809645d0 + h=0.5d0 + x1=x + y1=y + if (x.gt.3.5d0) x1=3.5d0 + if (y.gt.3.5d0) y1=3.5d0 + + xm=x1-y1 + xp=x1+y1 + xp2=-e+h*xp + zm=-Seq+z + +c Potential symmetrical in S1,S2,S3. Full set of +c derivatives for each order (1-4) and upper terms for others +C + v1=-3263.067522028298-9.17958228916955*xm**2- + y 0.5165611032611932*xm**4-0.5157212949525876*xm**6+ + y 18.31161578215203*xp2+ + y 30.14193751791963*xm**2*xp2-13.62543868575853*xm**4*xp2- + y 43.37232019119388*xp2**2-50.50353364079294*xm**2*xp2**2+ + y 70.36441193443143*xm**4*xp2**2+27.43935454999898*xp2**3+ + y 123.751990625258*xm**2*xp2**3-76.80240321256033*xm**4*xp2**3- + y 9.50017804016001*xp2**4-363.4487347625543*xm**2*xp2**4+ + y 113.1940248029751*xp2**5+376.6560011408163*xm**2*xp2**5- + y 164.6523756673548*xp2**6+9.16256842998227*zm + y -1.22230095639504*xm**2*zm+1.33032571356463*xp2* + y zm-0.94822119654751*xm**2*xp2*zm+0.7645470802285307*xm**4*xp2* + y zm-11.77270680473595*xp2**2*zm-0.4065994514809928*xm**4*xp2**2* + y zm-2.113651214829342*xp2**3*zm-3.653921741665064*xm**2*xp2**4* + y zm+26.53983199106825*xp2**6*zm+3.099164302936567*zm**2 + v2=-0.4668245990549825*xm**2*zm**2+0.05845413180128318*xm**4*zm**2 + y +2.708722250876111*xp2*zm**2-2.578482367020144*xm**2*xp2*zm**2- + y 0.1605392233404811*xm**4*xp2*zm**2-10.57780429022803*xp2**2* + y zm**2-3.496293826717189*xp2**3*zm**2+23.46280699747645*xp2**4* + y zm**2+1.8547816858377*zm**3-0.4003662844685243*xm**2*zm**3+ + y 3.040229985315839*xp2*zm**3-4.955739113923876*xp2**2* + y zm**3+14.05364889791468*xp2**3*zm**3-21.6926320924828*xp2**4* + y zm**3-1.321464834042384*zm**4+2.298844571392118*xm**2*zm**4- + y 2.633405421645483*xp2*zm**4+20.97178840867901*xp2**2* + y zm**4-32.18658937476802*xp2**3*zm**4-0.5992225949734171*zm**5+ + y 2.059827452250273*xp2*zm**5+0.6453850286056735*xp2**2* + y zm**5-0.4620689505336259*zm**6+0.7465042626807512*xp2* + y zm**6-0.1254018119377959*zm**7+0.01947721364782498*zm**8 + + VR=V1+V2 + +C SCALE AND SHIFT THE ZERO + VR=VR*1.0d-6 + RETURN + END + + SUBROUTINE PESleq6(V1,x,y,z) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + + V1=-3514.850376005703+1.189315215135138*(x-y)**2+ + - 0.824157459531989*(x-y)**4+0.03853108456851828*(x-y)**6+ + - 12.83265590340491*(-1.809659+0.5*(x+y))- + - 9.51736455454466*(x-y)**2*(-1.809659+0.5*(x+y))- + - 3.027576695974858*(x-y)**4*(-1.809659+0.5*(x+y))+ + - 10.94033338777717*(-1.809659+0.5*(x+y))**2+ + - 15.53332877554612*(x-y)**2*(-1.809659+0.5*(x+y))**2+ + - 6.063907309056958*(x-y)**4*(-1.809659+0.5*(x+y))**2- + - 13.79644533708051*(-1.809659+0.5*(x+y))**3- + - 26.67549601926293*(x-y)**2*(-1.809659+0.5*(x+y))**3- + - 6.200836894255189*(x-y)**4*(-1.809659+0.5*(x+y))**3+ + - 5.688103460541242*(-1.809659+0.5*(x+y))**4+ + - 52.9835500898771*(x-y)**2*(-1.809659+0.5*(x+y))**4- + - 11.88910471926647*(-1.809659+0.5*(x+y))**5- + - 43.99657824332825*(x-y)**2*(-1.809659+0.5*(x+y))**5+ + - 15.874512160015*(-1.809659+0.5*(x+y))**6- + - 8.60706112101134*(-1.824045+z)+ + - 1.264485336667462*(x-y)**2*(-1.824045+z)- + - 0.915127202947929*(-1.809659+0.5*(x+y))*(-1.824045+z) + v1=v1+ 0.6556566908758441*(x-y)**2*(-1.809659+0.5*(x+y))* + - (-1.824045+z)- + - 0.813078328219753*(x-y)**4*(-1.809659+0.5*(x+y))* + - (-1.824045+z)+ + - 12.42234678727481*(-1.809659+0.5*(x+y))**2* + - (-1.824045+z)+ + - 2.805488560712774*(x-y)**4*(-1.809659+0.5*(x+y))**2* + - (-1.824045+z)- + - 4.937250627623143*(-1.809659+0.5*(x+y))**3* + - (-1.824045+z)- + - 3.095201035295474*(x-y)**2*(-1.809659+0.5*(x+y))**4* + - (-1.824045+z)+ + - 18.85309150691318*(-1.809659+0.5*(x+y))**6* + - (-1.824045+z)- + - 3.209703208057476*(-1.824045+z)**2+ + - 0.5360421552708203*(x-y)**2*(-1.824045+z)**2- + - 0.263467844585989*(x-y)**4*(-1.824045+z)**2- + - 1.13298516075929*(-1.809659+0.5*(x+y))*(-1.824045+z)**2 + v1=v1- 0.06909229322445753*(x-y)**2*(-1.809659+0.5*(x+y))* + - (-1.824045+z)**2+ + - 1.649063526503709*(x-y)**4*(-1.809659+0.5*(x+y))* + - (-1.824045+z)**2+ + - 3.603611347474725*(-1.809659+0.5*(x+y))**2* + - (-1.824045+z)**2+ + - 3.757418764813337*(-1.809659+0.5*(x+y))**3* + - (-1.824045+z)**2+ + - 4.607672502246032*(-1.809659+0.5*(x+y))**4* + - (-1.824045+z)**2- + - 0.7490414640610651*(-1.824045+z)**3- + - 0.0888181500794012*(x-y)**2*(-1.824045+z)**3- + - 5.334303151299991*(-1.809659+0.5*(x+y))* + - (-1.824045+z)**3+ + - 1.37948603262339*(-1.809659+0.5*(x+y))**2* + - (-1.824045+z)**3+ + - 11.24395154910416*(-1.809659+0.5*(x+y))**3* + - (-1.824045+z)**3 - + - 17.85690001161674*(-1.809659+0.5*(x+y))**4*(-1.824045+z)**3 + v1=v1+ 0.7694433624551493*(-1.824045+z)**4- + - 0.939662303404418*(x-y)**2*(-1.824045+z)**4- + - 2.296000209594694*(-1.809659+0.5*(x+y))* + - (-1.824045+z)**4- + - 4.514249057965571*(-1.809659+0.5*(x+y))**2* + - (-1.824045+z)**4- + - 2.324765391545952*(-1.809659+0.5*(x+y))**3* + - (-1.824045+z)**4+ + - 0.223711667169141*(-1.824045+z)**5+ + - 1.164515013150094*(-1.809659+0.5*(x+y))* + - (-1.824045+z)**5- + - 2.825913168656484*(-1.809659+0.5*(x+y))**2* + - (-1.824045+z)**5+ + - 0.4811142779617512*(-1.824045+z)**6+ + - 1.292817090808966*(-1.809659+0.5*(x+y))* + - (-1.824045+z)**6+ + - 0.1657130839026308*(-1.824045+z)**7- + - 0.02192338698614548*(-1.824045+z)**8 + + v1=v1/1000000.0d0 + RETURN + END + + SUBROUTINE wifin(r1,r2,th,v) + +c +c In here almost all the points have been included, 342 in total +c +c + implicit real*8 (a-h,o-z) + + reoh=0.958649d0 + thetae=104.3475d0 + b1=2.0d0 + roh=0.951961d0 + alphaoh=2.587949757553683d0 + phh2=6.70164303995d0 + t0=0.01d0 + ut=20.d0 + x0=2.5d0 + ut2=20.d0 + x02=2.5d0 + + thetae=thetae*.314159265358979312d01*.00555555555555555555d0 + + + + + xs1=(r1+r2)*0.5d0-reoh + xs2=(r1-r2)*0.5d0 + xst=dcos(th)-dcos(thetae) + + + rs=dsqrt(0.5d0)*(r1+r2) + rm=dsqrt(0.5d0)*(r1-r2) + + rr1=r1-roh + rr2=r2-roh + + xep1=dexp(-2.d0*alphaoh*rr1)-2.d0*dexp(-alphaoh*rr1)+1.d0 + xep2=dexp(-2.d0*alphaoh*rr2)-2.d0*dexp(-alphaoh*rr2)+1.d0 + xep3=dexp(-b1*(rr1**2+rr2**2)) + rhh=dsqrt(r1**2+r2**2-2.d0*r1*r2*dcos(th)) + vhh=0.895714083584240987d6*dexp(-phh2*rhh) + vpb1=0.368181683997432046d5*xep1 + vpb2=0.368181683997432046d5*xep2 + + v0=-0.5710475632463d2 + vp1=-0.9552657956071d2*xst + y -0.6067577421565d4*xs1 + y +0.1857072188807d5*xst**2 + y -0.3804871940410d5*xs2**2 + y -0.2545869851885d5*xs1*xst + y -0.5068038862818d5*xs1**2 + y +0.1054898401242d4*xst**3 + y -0.2014577744857d5*xs2**2*xst + y -0.9931435344154d4*xs1*xst**2 + y +0.7065816455505d6*xs1*xs2**2 + y +0.2461185087692d5*xs1**2*xst + vp2=0.2211920749729d6*xs1**3 + y +0.3435869631455d4*xst**4 + y +0.5677448867330d5*xs2**2*xst**2 + y -0.4729349488906d6*xs2**4 + y -0.1062810010743d5*xs1*xst**3 + y -0.1711113900137d5*xs1*xs2**2*xst + y +0.6846686263643d5*xs1**2*xst**2 + y -0.2370095499969d7*xs1**2*xs2**2 + y -0.7435062393867d5*xs1**3*xst + y -0.5533282463806d6*xs1**4 + y -0.2653776328915d3*xst**5 + y -0.8787425252359d4*xs2**2*xst**3 + y -0.4792685061143d5*xs2**4*xst + vp3=0.7338375214709d4*xs1*xst**4 + y -0.2696329025096d5*xs1*xs2**2*xst**2 + y +0.3151302526351d7*xs1*xs2**4 + y -0.2680019540619d5*xs1**2*xst**3 + y -0.3411155916783d5*xs1**2*xs2**2*xst + y -0.2361101699812d5*xs1**3*xst**2 + y +0.4252318274899d7*xs1**3*xs2**2 + y +0.1909629119429d5*xs1**4*xst + y +0.1049781901349d7*xs1**5 + y +0.1007604315157d4*xst**6 + vp4=0.3955134247500d5*xs2**2*xst**4 + y -0.1254283199622d7*xs2**6 + y -0.4742047112461d4*xs1*xst**5 + y -0.1237633931215d6*xs1*xs2**2*xst**3 + y +0.1989124023766d5*xs1**2*xst**4 + y +0.3708558366230d6*xs1**2*xs2**2*xst**2 + y +0.6018342068259d5*xs1**4*xst**2 + y -0.1239376779228d7*xs1**6 + y +0.8543568929153d4*xs2**2*xst**5 + y -0.8470361074768d5*xs2**6*xst + y -0.4323120218548d4*xs1*xst**6 + y -0.7691197466131d5*xs1*xs2**2*xst**4 + y -0.4892966870779d7*xs1*xs2**6 + vp5=0.4364087781854d5*xs1**2*xst**5 + y -0.3909155940947d8*xs1**3*xs2**4 + y -0.1712950976860d7*xs1**4*xs2**2*xst + y +0.2087371244764d6*xs1**5*xst**2 + y -0.2006710729271d8*xs1**5*xs2**2 + y +0.8789956427347d6*xs1**7 + y -0.4317629382959d5*xs2**2*xst**6 + y +0.1618789634550d6*xs2**4*xst**4 + y +0.1587687527646d7*xs2**6*xst**2 + y +0.3990003280266d7*xs2**8 + y +0.2817080082968d6*xs1*xs2**2*xst**5 + y +0.4482768991684d6*xs1*xs2**4*xst**3 + y -0.4070569279958d5*xs1**2*xst**6 + vp6=-0.8864550152825d6*xs1**2*xs2**2*xst**4 + y +0.6290955366349d8*xs1**2*xs2**6 + y +0.8993645659488d5*xs1**3*xst**5 + y +0.9210094589981d6*xs1**3*xs2**2*xst**3 + y -0.1753794544813d6*xs1**4*xst**4 + y +0.1311789234505d9*xs1**4*xs2**4 + y +0.1060664061203d6*xs1**5*xst**3 + y +0.5340569244466d7*xs1**5*xs2**2*xst + y -0.3105624184638d6*xs1**6*xst**2 + vp7=0.4571723553398d8*xs1**6*xs2**2 + y +0.7470215166875d6*xs1**7*xst + y -0.2218921351107d6*xs2**4*xst**5 + y +0.1316232741184d5*xs1*xst**8 + y +0.2023964877189d6*xs1*xs2**2*xst**6 + y -0.7539991140488d7*xs1*xs2**6*xst**2 + y -0.3438106580968d8*xs1*xs2**8 + y -0.1065626654616d6*xs1**2*xst**7 + y -0.2295235891775d6*xs1**2*xs2**2*xst**5 + y +0.1796585516193d6*xs1**3*xst**6 + y -0.2005145616345d7*xs1**3*xs2**4*xst**2 + vp8=-0.1375861482407d9*xs1**3*xs2**6 + y -0.3423607622386d6*xs1**4*xst**5 + y +0.6119673444298d7*xs1**4*xs2**2*xst**3 + y -0.1085580079570d8*xs1**5*xs2**2*xst**2 + y -0.1807866589974d9*xs1**5*xs2**4 + y -0.4710714407573d7*xs1**6*xs2**2*xst + y -0.3555338273642d8*xs1**7*xs2**2 + y -0.4225310291748d7*xs1**8*xst + y +0.6802694583683d5*xs2**2*xst**8 + y -0.3612806511976d6*xs2**6*xst**4 + y -0.1984202384497d7*xs2**8*xst**2 + vp9=-0.5762585608683d4*xs1*xst**9 + y -0.4190242262127d6*xs1*xs2**2*xst**7 + y -0.6862320609996d6*xs1*xs2**4*xst**5 + y +0.6897610937412d5*xs1**2*xst**8 + y +0.1979393850517d7*xs1**2*xs2**2*xst**6 + y +0.1886232859991d7*xs1**2*xs2**4*xst**4 + y +0.1704308445329d8*xs1**2*xs2**6*xst**2 + y +0.5734113225539d8*xs1**2*xs2**8 + y -0.1909984371034d6*xs1**3*xst**7 + y -0.4274427341543d7*xs1**3*xs2**2*xst**5 + y -0.3010333878947d7*xs1**3*xs2**4*xst**3 + y +0.2957770502057d6*xs1**4*xst**6 + vp10=0.4275754424910d7*xs1**4*xs2**2*xst**4 + y +0.7954251656913d8*xs1**4*xs2**6 + y -0.1018371344211d8*xs1**5*xs2**2*xst**3 + y +0.2825744251920d6*xs1**6*xst**4 + y +0.1682459960715d8*xs1**6*xs2**2*xst**2 + y +0.1012584539651d9*xs1**6*xs2**4 + y +0.5160076143117d7*xs1**9*xst + y -0.4748852759269d3*xst**11 + y -0.1946436003678d5*xs2**2*xst**9 + + vp=vp1+vp2+vp3+vp4+vp5+vp6+vp7+vp8+vp10+vp9 + + vps1=42395.535333d0*xep1 + vps2=42395.535333d0*xep2 + + y1=1.d0/(1.d0+dexp(ut*(x0-r1))) + y2=1.d0/(1.d0+dexp(ut*(x0-r2))) + y12=1.d0/(1.d0+dexp(ut2*(x02-r1))) + y22=1.d0/(1.d0+dexp(ut2*(x02-r2))) + + vp=vp*xep3*(1-y12)*(1-y22) + voh1=vpb1*(1-y1)+y1*vps1 + voh2=vpb2*(1-y2)+y2*vps2 + + v=v0+vp+voh1+voh2+vhh + + return + end + + + SUBROUTINE cvps(V,rij1,rij2,rij3) +c subroutine vibpot(rij,v,n) + implicit real*8 (a-h,o-z) +c +c pes for h2o, +c Harry Partridge and David W. Schwenke, J. Chem. Phys., +c submitted Nov. 8, 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 mass dependent factors are included. the nuclear masses +c should be passed to this program using the array xm in +c common potmcm. xm(1) is the +c mass of the hydrogen associated with rij(i,1), and xm(2) +c is the mass of the hydrogen associated with rij(i,2). +c all masses are in au. +c + integer :: i, j, idx, idxm, ifirst + dimension c5z(245),cbasis(245),ccore(245), + $ crest(245),idx(245,3),fmat(15,3),cmass(9),idxm(9,3) +c common/potrot/fact1,fact2,c1,s1,icoord,xm(2),xmx,iperm +c $$$ common/potmcm/xm(2) +c +c expansion indicies +c + 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/ +c +c expansion coefficients for 5z ab initio data +c + data (c5z(i),i=1,45)/ + $ 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/ + + data (c5z(i),i=46,90)/ + $-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/ + + data (c5z(i),i=91,135)/ + $-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/ + + + data (c5z(i),i=136,180)/ + $ 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/ + + data (c5z(i),i=181,225)/ + $ 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/ + + + data (c5z(i),i=226,245)/ + $ 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/ +c +c expansion coefficients for basis correction +c + data (cbasis(i),i=1,45)/ + $ 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/ + + + data (cbasis(i),i=46,90)/ + $-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/ + + + + data (cbasis(i),i=91,135)/ + $ 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/ + + + + + data (cbasis(i),i=136,180)/ + $ 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/ + + + data (cbasis(i),i=181,225)/ + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 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 (cbasis(i),i=226,245)/ + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00/ +c +c expansion coefficients for core correction +c + data (ccore(i),i=1,45)/ + $ 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/ + + + data (ccore(i),i=46,90)/ + $-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/ + + + data (ccore(i),i=91,135)/ + $-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/ + + + data (ccore(i),i=136,180)/ + $ 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/ + + data (ccore(i),i=181,225)/ + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 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=226,245)/ + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00/ +c +c expansion coefficients for v rest +c + data (crest(i),i=1,45)/ + $ 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/ + + data (crest(i),i=46,90)/ + $ 0.0000000000000D+00, 0.0000000000000D+00, 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/ + + data (crest(i),i=91,135)/ + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 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=136,180)/ + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 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=181,225)/ + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 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=226,245)/ + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, + $ 0.0000000000000D+00, 0.0000000000000D+00/ +c +c expansion indicies for mass correction +c + data idxm/1,2,1,1,3,2,1,2,1, + $ 2,1,1,3,1,2,2,1,1, + $ 1,1,2,1,1,1,2,2,3/ +c +c expansion coefficients for mass correction +c + data cmass/ -8.3554183D+00,3.7036552D+01,-5.2722136D+00, + $ 1.6843857D+01,-7.0929741D+01,5.5380337D+00,-2.9962997D+01, + $ 1.3637682D+02,-3.0530195d+00/ +c +c two body parameters +c + 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/ +c +c scaling factors for contributions to emperical potential +c +c data f5z,fbasis,fcore,frest/0.99967788500000d0, +c $ 0.15860145369897d0,-1.6351695982132d0,1d0/ + data f5z,fbasis,fcore,frest/0.0d0,0.0d0,-1.0d0,0.0d0/ + save + data ifirst/0/ + if(ifirst.eq.0)then + ifirst=1 +c write(6,1) + 1 format(/1x,'pes for h2o', + $ /1x,'by Harry Partridge and David W. Schwenke', + $ /1x,'submitted to J. Chem. Phys. Nov. 8, 1996') +c write(6,56) + 56 format(/1x,'parameters before adjustment') +c 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) +c 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) + 2 continue +c +c remove mass correction from vrest +c + xmh=1836.152697d0 + xmhi=1d0/xmh + xmd=3670.483031d0 + fact=1.d0/((1.d0/xmd)-(1.d0/xmh)) +c write(6,65) + 65 format(/1x,'parameters for delta v hdo ', + $ /1x,' i j k') + do 60 i=1,9 +c write(6,5)(idxm(i,j)-1,j=1,3),cmass(i) + cmass(i)=cmass(i)*fact + corr=cmass(i)*xmhi + if(idxm(i,1).eq.idxm(i,2))corr=corr*0.5d0 + do 61 j=1,245 + if(idx(j,1).eq.idxm(i,1).and.idx(j,2).eq.idxm(i,2).and. + $ idx(j,3).eq.idxm(i,3))then + crest(j)=crest(j)-corr + go to 62 + end if + 61 continue + 62 continue + do 63 j=1,245 + if(idx(j,2).eq.idxm(i,1).and.idx(j,1).eq.idxm(i,2).and. + $ idx(j,3).eq.idxm(i,3))then + crest(j)=crest(j)-corr + go to 64 + end if + 63 continue + 64 continue + 60 continue +c write(6,70)xm +c 70 format(/1x,'masses used for mass correction: ',1p2e15.7) + xm1=1.d0/1836.152697d0 + xm2=1.d0/1836.152697d0 +c +c adjust parameters using scale factors +c +c 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 + do 59 i=1,245 + c5z(i)=f5z*c5z(i)+fbasis*cbasis(i)+fcore*ccore(i) + $ +frest*crest(i) + 59 continue +c 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) + 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)) + do 66 i=1,9 + cmass(i)=cmass(i)*frest + 66 continue +c write(6,76)((idxm(i,j),j=1,3),cmass(i),i=1,9) + 76 format(/1x,'mass correction factors ', + $ /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 + do 67 i=1,9 + cmass(i)=cmass(i)*4.556335d-6 + 67 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 +c do 6 i=1,n + x1=(rij1-reoh)/reoh + x2=(rij2-reoh)/reoh + x3=cos(rij3)-ce + rhh=sqrt(rij1**2+rij2**2 + $ -2d0*rij1*rij2*cos(rij3)) + vhh=phh1*exp(-phh2*rhh) + ex=exp(-alphaoh*(rij1-roh)) + voh1=deoh*ex*(ex-2d0) + ex=exp(-alphaoh*(rij2-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=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=v+term + 12 continue + v1=0d0 + v2=0d0 + do 13 j=1,9 + v1=v1+cmass(j)*fmat(idxm(j,1),1)*fmat(idxm(j,2),2) + $ *fmat(idxm(j,3),3) + v2=v2+cmass(j)*fmat(idxm(j,2),1)*fmat(idxm(j,1),2) + $ *fmat(idxm(j,3),3) + 13 continue + v=v+xm1*v1+xm2*v2 + +c write(6,*)v,xm1,v1,xm2,v2 + + v=v*exp(-b1*((rij1-reoh)**2+(rij2-reoh)**2)) + $ +c5z(1) + $ +voh1+voh2+vhh +c 6 continue + + + return + end diff --git a/tests/H2O_CVRQD/ERROR.ref b/tests/H2O_CVRQD/ERROR.ref new file mode 100644 index 00000000..783ca168 --- /dev/null +++ b/tests/H2O_CVRQD/ERROR.ref @@ -0,0 +1 @@ +ERROR in force_h2o.F90: Numerical forces not yet implemented! diff --git a/tests/H2O_CVRQD/input.in b/tests/H2O_CVRQD/input.in new file mode 100644 index 00000000..aacf36fb --- /dev/null +++ b/tests/H2O_CVRQD/input.in @@ -0,0 +1,19 @@ +! Test the analytical CVRQD water potential + +&general +nstep=1, +irest=0, +idebug=1 + +pot='_h2o_' +h2opot='cvrqd' +mdtype='MD', ! classical MD +dt=20., ! number of steps and timestep +nstep=2 +/ + +&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_CVRQD/mini.xyz b/tests/H2O_CVRQD/mini.xyz new file mode 100644 index 00000000..85d61ae3 --- /dev/null +++ b/tests/H2O_CVRQD/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 diff --git a/tests/test.sh b/tests/test.sh index b4f57a0b..97a2d96e 100755 --- a/tests/test.sh +++ b/tests/test.sh @@ -120,6 +120,8 @@ restart_sh.bin restart_sh.bin.old restart_sh.bin.?? restart.xyz.old restart.xyz. # Run all tests if [[ $TESTS = "all" ]];then + # TODO: Automaticall select all folders except for the special cases + # that depend on optional features. folders=(INIT CMD NHC-GLOBAL SHAKE \ SH_EULER SH_RK4 SH_BUTCHER SH_RK4_PHASE \ SH_IGNORE SH_NACM_FAIL SH_S0S1 SH_ENERGY_DIFF SH_ENERGY_DRIFT \ @@ -127,7 +129,7 @@ if [[ $TESTS = "all" ]];then LZ_SS LZ_ST LZ_ENE \ PIMD ABINITIO ABINITIO-FAIL MTS \ LANGEVIN QT QT2 PIGLE PIGLE2 GLE-CANONICAL \ - H2O_SCHWENKE HARMON MORSE DOUBLEWELL SPLINE MM MINI QMMM \ + H2O_SCHWENKE H2O_CVRQD HARMON MORSE DOUBLEWELL SPLINE MM MINI QMMM \ ANALYZE_EXT CMDLINE WATER_FAIL ERMD) if [[ $MPI = "TRUE" ]];then