diff --git a/ocaml/Input.ml b/ocaml/Input.ml index 6a1ff43..70b321b 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -817,13 +817,6 @@ end = struct let _ = Lazy.force Qputils.ezfio_filename in - let () = - match (Pseudo.read () |> Pseudo.to_bool, t) with - | (false, _) - | (true , None) -> () - | _ -> failwith "Jastrow and Pseudopotentials are incompatible for now" - in - to_string t |> Ezfio.set_jastrow_jast_type @@ -902,8 +895,6 @@ let validate () = Sampling.read () and ts = Time_step.read () - and jast_type = - Jastrow_type.read () and do_pseudo = Pseudo.read () in @@ -955,14 +946,6 @@ let validate () = | _ -> () in - (* Pseudo and Jastrow are incompatible *) - let () = - match (Pseudo.to_bool do_pseudo, jast_type) with - | (true, Jastrow_type.Core ) - | (true, Jastrow_type.Simple) -> failwith "Jastrow and Pseudopotentials are incompatible" - | _ -> () - in - (* Fitcusp is incompatible with pseudo *) let () = let f = diff --git a/src/pseudo.irp.f b/src/pseudo.irp.f index dd01e01..156c25a 100644 --- a/src/pseudo.irp.f +++ b/src/pseudo.irp.f @@ -333,31 +333,54 @@ real function ylm(l,m,x,y,z,r_inv) continue case(1) + ylm = ylm * sq3 * r_inv select case (m) case (-1) - ylm = ylm * sq3 * y * r_inv + ylm = ylm * y case (0) - ylm = ylm * sq3 * z * r_inv + ylm = ylm * z case (1) - ylm = ylm * sq3 * x * r_inv + ylm = ylm * x + end select + + case(2) + ylm = ylm * r_inv * r_inv * sqrt(5.) + select case (m) + case(-2) + ylm = ylm * sqrt(3.) * x * y + case(-1) + ylm = ylm * sqrt(3.) * y * z + case(0) + ylm = ylm * 0.5 * (2.*z*z - x*x - y*y) + case(1) + ylm = ylm * sqrt(3.) * z * x + case(2) + ylm = ylm * 0.5 * sqrt(3.) * (x*x - y*y) + end select + + case(3) + ylm = ylm * r_inv * r_inv * r_inv * sqrt(7.) + select case (m) + case(-3) + ylm = ylm * 0.25 * sqrt(10.) * (3.*x*x - y*y) + case(-2) + ylm = ylm * sqrt(15.) * x*y*z + case(-1) + ylm = ylm * 0.25*sqrt(6.) * y * (4.*z*z - x*x -y*y) + case(0) + ylm = ylm * 0.5 * z * (2.*z*z - 3.*(x*x + y*y)) + case(1) + ylm = ylm * 0.25*sqrt(6.) * x * (4.*z*z - x*x -y*y) + case(2) + ylm = ylm * 0.5*sqrt(15.) * z * (x*x -y*y) + case(3) + ylm = ylm * 0.25*sqrt(10.) * x * (x*x - 3. * y*y) end select -! case(2) -! select case (m) -! case(-2) -! ylm = ylm * sqrt(15.) * x * y * r_inv * r_inv -! case(-1) -! ylm = ylm * sqrt(15.) * y * z * r_inv * r_inv -! case(0) -! ylm = 0.5 * ylm * sqrt(15.) * (2.*z*z - x*x - y*y) * r_inv * r_inv -! case(1) -! ylm = ylm * sqrt(15.) * z * x * r_inv * r_inv -! case(2) -! ylm = 0.5 * ylm * sqrt(15.) * (x*x - y*y) * r_inv * r_inv -! end select case default - stop 'problem in Ylm of pseudo' + print *, 'l=', l + stop 'problem in Ylm of pseudo : Ylm not implemented (pseudo.irp.f)' end select