Working code version (#1)

.github/workflows/CI.yml

@@ -5,6 +5,8 @@ on:
       - main
     tags: ['*']
   pull_request:
+env:
+  JULIA_NUM_THREADS: 2
 concurrency:
   # Skip intermediate builds: always.
   # Cancel intermediate builds: only if it is a pull request build.
@@ -31,5 +33,9 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v1
+      - name: MatrixLabTools repo install
+        run: julia --color=yes --check-bounds=yes -e 'using Pkg;
+               Pkg.Registry.add(RegistrySpec(url="https://github.com/MatrixLabTools/PackageRegistry"));
+               Pkg.Registry.add(RegistrySpec(url="https://github.com/JuliaRegistries/General"));'
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1

.gitignore

@@ -1 +1,4 @@
-/Manifest.toml
+Manifest.toml
+.CondaPkg
+*.clean
+*.dat

CondaPkg.toml

@@ -0,0 +1,5 @@
+channels = ["anaconda", "conda-forge", "psi4", "conda-forge/label/libint_dev"]
+
+[deps]
+psi4 = ">=1.4"
+pytest = ""

Project.toml

@@ -3,8 +3,14 @@ uuid = "e2b7ad92-acfe-472f-ae48-51c2aed5b591"
 authors = ["Teemu Järvinen <teemu.j.jarvinen@gmail.com> and contributors"]
 version = "0.1.0-DEV"
 
+[deps]
+AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
+PotentialCalculation = "76415700-ec45-11e8-0f65-f1df3c899b21"
+PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+
 [compat]
 julia = "1.9"
+PotentialCalculation = "0.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -1,3 +1,31 @@
 # Psi4Calculator
 
 [![Build Status](https://github.com/MatrixLabTools/Psi4Calculator.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/MatrixLabTools/Psi4Calculator.jl/actions/workflows/CI.yml?query=branch%3Amain)
+
+This is a [Psi4](https://psicode.org) calculator for [PotentialCalculation](https://github.com/MatrixLabTools/PotentialCalculation.jl)
+
+## Installation
+
+Hit "]" to enter "pkg>"
+
+```julia
+pkg> add registry add https://github.com/MatrixLabTools/PackageRegistry
+pkg> add https://github.com/MatrixLabTools/Psi4Calculator.jl
+```
+
+## Example use case
+
+```julia
+using PotentialCalculation
+using Psi4Calculator
+
+N2 = isolated_system( [Atom(:N, [1., 0., 0.].*u"Å"), Atom(:N, [0., 0., 0.].*u"Å")] )
+
+ca = Calculator("blyp", "def2-svp",Psi4(memory="1000MiB", nthreads=2))
+
+calculate_energy(ca, N2)
+```
+
+## Note
+
+This is still experimental and may or may not work.

src/Psi4Calculator.jl

@@ -1,5 +1,91 @@
 module Psi4Calculator
 
 # Write your package code here.
+using AtomsBase
+using PythonCall
+using PotentialCalculation.Calculators
+using PotentialCalculation.Clusters
+
+export Psi4
+
+"""
+    gpsi4init=false
+
+Global variable to see if Psi4 was initiated for current process
+"""
+global gpsi4init=false
+
+"""
+    gPsi4 = undef
+
+Hold Psi4 object given by PyCall. Used to do calculation in current process
+"""
+global gPsi4 = undef
+
+"""
+    initpsi4(;memory="500 MiB", quiet=true, nthreads=1)
+
+Used to intialize Psi4 environment
+"""
+function initpsi4(;memory="500 MiB", quiet=true, nthreads=1)
+    global gPsi4 = pyimport("psi4")
+    global gpsi4init = true
+    quiet && gPsi4.core.be_quiet()
+    gPsi4.set_memory(memory)
+    nthreads > 1 && gPsi4.set_num_threads(nthreads)
+end
+
+"""
+    mutable struct Psi4 <: AbstractCalculationProgram
+
+Holds information that calculations are to be done with Psi4
+
+# Fields
+- `memory="500MiB"` : memory used by Psi4
+- `nthreads=1`  : number of threads used by Psi4
+"""
+mutable struct Psi4 <: AbstractCalculationProgram
+    memory::String
+    nthreads::UInt
+    function Psi4(;memory="500MiB", nthreads=1)
+        initpsi4(memory=memory, nthreads=nthreads)
+        @debug gPsi4
+        new(memory, nthreads)
+    end
+end
+
+
+function Calculators.calculate_energy(cal::Calculator{Psi4}, point::Union{Cluster,AbstractSystem};
+                                      basename="base", ghost=undef, id="", pchannel=undef)
+    ! gpsi4init && initpsi4(memory=cal.calculator.memory, nthreads=cal.calculator.nthreads)
+    s=sprint( (io, x) -> print_xyz(io,x, printheader=false), point)
+    c = gPsi4.geometry(s)
+    out = gPsi4.energy(cal.method*"/"*cal.basis, molecule=c)
+    pchannel != undef && put!(pchannel,true)
+    return pyconvert(Float64, out)
+end
+
+
+function Calculators.calculate_energy(cal::Calculator{Psi4}, points;
+                                      basename="base", ghost=undef, id="", pchannel=undef)
+    return map( x -> calculate_energy(cal, x, basename=basename, ghost=ghost, id=id, pchannel=pchannel), points )
+end
+
+
+function Calculators.bsse_corrected_energy(cal::Calculator{Psi4}, c1::Union{Cluster,AbstractSystem}, c2::Union{Cluster,AbstractSystem};
+                                           basename="base", id="", pchannel=undef)
+    ! gpsi4init && initpsi4(memory=cal.calculator.memory)
+    s1=sprint( (io, x) -> print_xyz(io,x, printheader=false), c1)
+    s2=sprint( (io, x) -> print_xyz(io,x, printheader=false), c2)
+    c = gPsi4.geometry(s1*"\n--\n"*s2)
+    out = gPsi4.energy(cal.method*"/"*cal.basis, molecule=c, bsse_type="cp")
+    pchannel != undef && put!(pchannel,true)
+    return pyconvert(Float64, out)
+end
+
+function Calculators.bsse_corrected_energy(cal::Calculator{Psi4}, c1, c2;
+                                           basename="base", id="", pchannel=undef)
+    return map( (x,y) -> bsse_corrected_energy(cal, x, y, basename=basename, id=id, pchannel=pchannel ), c1, c2  )
+end
 
 end

test/Project.toml

@@ -0,0 +1,4 @@
+[deps]
+Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+PotentialCalculation = "76415700-ec45-11e8-0f65-f1df3c899b21"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/runtests.jl

@@ -1,6 +1,49 @@
-using Psi4Calculator
 using Test
 
+using Distributed
+addprocs(2)
+
+@everywhere using PotentialCalculation
+@everywhere using Psi4Calculator
+
+fname = tempname() * ".jld2"
+rname = tempname() * ".jld2"
+sname = tempname() * ".jld2"
+xyzname = tempname() * ".xyz"
+
+formic_acid=Cluster(
+ [-6.7041359778      1.3501192944      0.0102209137
+-5.3688853815      1.2229556023      0.0440598937
+-7.2470157373      2.4374213225      0.0651311769
+ -5.0398812618      2.1435406993      0.1155201154
+ -7.2001330967      0.3718768293     -0.0703451879]',
+  AtomOnlySymbol.(["C", "O", "O", "H", "H"]) )
+
+Ar = Cluster( Atom(:Ar, rand(3)u"Å"  )  )
+N2 = isolated_system( [Atom(:N, [1., 0., 0.].*u"Å"), Atom(:N, [0., 0., 0.].*u"Å")] )
+
+open(xyzname,"w") do io
+    print_xyz(io, formic_acid)
+end
+
+pbar=true
+
+testrestarts = false
+
 @testset "Psi4Calculator.jl" begin
-    # Write your tests here.
+    ca = Calculator("blyp", "def2-svp",Psi4(memory="1000MiB", nthreads=2))
+
+    input1=create_inputs(xyzname, Ar, ca)
+    inputs=create_inputs(xyzname, N2, ca; npoints=5)
+    inputss=create_inputs(xyzname, xyzname, ca)
+
+    data1=calculate_potential(inputs, save_file=fname, pbar=pbar)
+    data2=calculate_potential(fname,ca,save_file=sname, restart_file=rname, pbar=pbar)
+    data3=continue_calculation(rname,ca, save_file=sname, restart_file=rname, pbar=pbar)
+
+    calculate_energy(ca, N2)
+    calculate_energy(ca, [N2,N2])
+    @test all(isapprox.(data1["Energy"], data2["Energy"], atol=2E-6))
+    @test all(isapprox.(data1["Energy"], data3["Energy"], atol=2E-6))
+    @test all(isapprox.(data2["Energy"], data3["Energy"], atol=2E-6))
 end