Merge pull request #1 from fkfest/auto/from-devel-repo

Merge changes from development repository

.github/workflows/push-to-stable.yml

@@ -0,0 +1,63 @@
+# @author jphaupt
+name: Push changes to a public repository
+
+on:
+  push:
+    branches:
+      - stable # branch on private that triggers action
+
+# edit these environment variables if you wish to use in a new project
+# note this depends on a separate workflow called build_and_test to exist
+# (this job does not run unless build_and_test is successful)
+env:
+  BOT_TOKEN: ${{ secrets.PRIVATE_TO_PUBLIC_PAT }} # PAT for jph-bot
+  TARGET_REPO: e-cojl
+  TARGET_BRANCH: main
+  SOURCE_REPO: e-cojl-devel
+  SOURCE_BRANCH: stable
+
+jobs:
+  push-changes:
+    runs-on: ubuntu-latest
+    steps:
+      - name: checkout code
+        uses: actions/checkout@v2
+        with:
+          ref: stable # unfortunately, cannot use env here, must change manually as well
+          fetch-depth: 0
+
+      - name: clone $SOURCE_BRANCH repo
+        run: |
+          git clone https://fkfest:[email protected]/fkfest/$TARGET_REPO.git
+          cd $TARGET_REPO
+          git checkout -b auto/from-devel-repo
+          git config user.name "actions-user"
+          git config user.email "[email protected]"
+      - name: pull changes into local copy of public repo from private repo
+        run: |
+          cd $TARGET_REPO
+          git remote add source https://fkfest:[email protected]/fkfest/$SOURCE_REPO
+          git fetch source $SOURCE_BRANCH
+      - name: merge changes in public repo
+        run: |
+          cd $TARGET_REPO
+          git merge source/$SOURCE_BRANCH --no-edit
+      - name: publish new branch
+        run: |
+          cd $TARGET_REPO
+          git push -u origin auto/from-devel-repo
+      - name: Create Pull Request
+        run: |
+          cd $TARGET_REPO
+          PR_TITLE="Merge changes from development repository"
+          PR_BODY="This pull request merges changes from the private development repository and was automatically generated by GitHub Actions."
+          PR_HEAD="auto/from-devel-repo" # branch with changes
+          PR_BASE="$TARGET_BRANCH"
+          PR_API_URL=$(curl --silent --request POST --url https://api.github.com/repos/fkfest/$TARGET_REPO/pulls \
+            --header "authorization: Bearer $BOT_TOKEN" \
+            --header "content-type: application/json" \
+            --data "{\"title\":\"$PR_TITLE\",\"body\":\"$PR_BODY\",\"head\":\"$PR_HEAD\",\"base\":\"$PR_BASE\"}" \
+            --write-out '%{url_effective}' \
+            --output /dev/null)
+          PR_NUMBER=$(echo $PR_API_URL | cut -d / -f 8)
+          echo "Pull Request #$PR_NUMBER created."

README.md

@@ -26,3 +26,22 @@ Various options are available (use `-h` option for a list of `e-co.jl` options):
 
 Default scratch dir is `./ecojlscr` and default fcidump file is `FCIDUMP`.
 
+```
+Electron coil
+A poem by Bing
+
+In the heart of an atom lies a tiny core
+Where protons and neutrons are tightly bound
+But around this nucleus, there's so much more
+A cloud of electrons that swirls around
+
+They don't orbit in circles like planets do
+But jump and spin in quantum states
+They can be here and there and everywhere too
+And sometimes they even change their mates
+
+This is the electron coil, the turmoil of the shell
+The source of light and heat and power
+The force that makes the atoms repel or gel
+The spark that ignites the cosmic flower
+```

src/dfhf.jl

@@ -1,12 +1,12 @@
 module DFHF
 using LinearAlgebra, TensorOperations, Printf
 using ..ECInfos
-using ..MSystem
 using ..ECInts
+using ..MSystem
 using ..DIIS
 using ..TensorTools
 
-export dfhf
+export dfhf, GuessType, GUESS_HCORE, GUESS_SAD
 
 """ integral direct df-hf """ 
 function dffock(EC,cMO,bao,bfit)
@@ -103,20 +103,63 @@ function generate_integrals(ms::MSys, EC::ECInfo; save3idx = true)
   else
     save(EC,"PL",M)
   end
-  return ECInts.nuclear_repulsion(bao)
+  return nuclear_repulsion(ms)
+end
+
+@enum GuessType GUESS_HCORE GUESS_SAD GUESS_GWH GUESS_ORB
+
+function guess_hcore(EC::ECInfo)
+  hsmall = load(EC,"hsmall")
+  sao = load(EC,"sao")
+  ϵ,cMO = eigen(Hermitian(hsmall),Hermitian(sao))
+  return cMO
+end
+
+function guess_sad(ms::MSys, EC::ECInfo)
+  # minao = "ano-rcc-mb"
+  minao = "ano-r0"
+  # minao = "sto-6g"
+  bminao = BasisSet(minao,genxyz(ms,bohr=false))
+  bao,bfit = generate_basis(ms)
+  smin2ao = overlap(bminao,bao)
+  eldist = electron_distribution(ms,minao)
+  saoinv = invchol(Hermitian(load(EC,"sao")))
+  # display(eldist)
+  denao = saoinv * smin2ao' * diagm(eldist) * smin2ao * saoinv
+  # dc = nc
+  n,cMO = eigen(Hermitian(-denao))
+  # display(n)
+  return cMO
+end
+
+function guess_gwh(ms::MSys, EC::ECInfo)
+  error("not implemented yet")
+end
+
+function guess_orb(ms::MSys, EC::ECInfo, guess::GuessType)
+  if guess == GUESS_HCORE
+    return guess_hcore(EC)
+  elseif guess == GUESS_SAD
+    return guess_sad(ms,EC)
+  elseif guess == GUESS_GWH
+    return guess_gwh(ms,EC)
+  elseif guess == GUESS_ORB
+    return load(EC,"cMO")
+  else
+    error("unknown guess type")
+  end
 end
 
-function dfhf(ms::MSys, EC::ECInfo; direct = false)
+function dfhf(ms::MSys, EC::ECInfo; direct = false, guess = GUESS_SAD)
   diis = Diis(EC.scr)
   thren = sqrt(EC.thr)*0.1
   Enuc = generate_integrals(ms, EC; save3idx=!direct)
   if direct
     bao,bfit = generate_basis(ms)
   end
+  cMO = guess_orb(ms,EC,guess)
   hsmall = load(EC,"hsmall")
   sao = load(EC,"sao")
-  ϵ,cMO = eigen(Hermitian(hsmall),Hermitian(sao))
-  # display(ϵ)
   SP = EC.space
   previousEHF = 0.0
   println("Iter     Energy      DE          Res         Time")

src/dump.jl

@@ -80,7 +80,7 @@ function read_fcidump(fcidump::String)
   if simtra
     println("Non-Hermitian")
   end
-  if isnothing(headvar(fd, "NPY2")) && isnothing(headvar(fd, "NPYAA"))
+  if isnothing(headvar(fd, "NPY2")) && isnothing(headvar(fd, "NPY2AA"))
     # read integrals from fcidump file
     read_integrals!(fd,fdf)
     close(fdf)
@@ -105,7 +105,7 @@ function read_header(fdfile)
     if length(line) == 0
       continue
     end
-    if line == "/"
+    if line == "/" || line == "&END"
       # end of header
       break
     end
@@ -127,7 +127,8 @@ function read_header(fdfile)
             # store the previous array
             head[variable_name] = elements
           end
-          variable_name = prev_el
+          # case-insensitive variable names in the header
+          variable_name = uppercase(prev_el)
           elements = []
           prev_el = ""
         else
@@ -164,11 +165,11 @@ function read_integrals!(fd::FDump, dir::AbstractString)
     fd.int2 = mmap_integrals(fd, dir, "NPY2")
     fd.int1 = mmap_integrals(fd, dir, "NPY1")
   else
-    fd.int2aa = mmap_integrals(fd, dir, "NPYAA")
-    fd.int2bb = mmap_integrals(fd, dir, "NPYBB")
-    fd.int2ab = mmap_integrals(fd, dir, "NPYAB")
-    fd.int1a = mmap_integrals(fd, dir, "NPYA")
-    fd.int1b = mmap_integrals(fd, dir, "NPYB")
+    fd.int2aa = mmap_integrals(fd, dir, "NPY2AA")
+    fd.int2bb = mmap_integrals(fd, dir, "NPY2BB")
+    fd.int2ab = mmap_integrals(fd, dir, "NPY2AB")
+    fd.int1a = mmap_integrals(fd, dir, "NPY1A")
+    fd.int1b = mmap_integrals(fd, dir, "NPY1B")
   end
   if isnothing(headvar(fd, "ENUC"))
     error("ENUC option not found in fcidump")
@@ -280,7 +281,8 @@ function read_integrals!(fd::FDump, fdfile::IOStream)
     line = split(linestr)
     if length(line) != 5
       # println("Last line: ",linestr)
-      break
+      # skip lines (in the case there is something left from header)...
+      continue
     end
     integ = parse(Float64,line[1])
     i1 = parse(Int,line[2])
@@ -356,4 +358,4 @@ function mmap_integrals(fd::FDump, dir::AbstractString, key::AbstractString)
   return mnpymmap(file)
 end
 
-end #module
+end #module

src/e-co.jl

@@ -102,10 +102,11 @@ function main()
   println(size(EC.fd.int2))
   norb = headvar(EC.fd, "NORB")
   nelec = headvar(EC.fd, "NELEC")
+  ms2 = headvar(EC.fd, "MS2")
 
   SP = EC.space
 
-  SP['o'], SP['v'], SP['O'], SP['V'] = get_occvirt(EC, occa, occb, norb, nelec)
+  SP['o'], SP['v'], SP['O'], SP['V'] = get_occvirt(EC, occa, occb, norb, nelec, ms2)
   SP[':'] = 1:headvar(EC.fd,"NORB")
 
 

src/ecinfos.jl

@@ -70,7 +70,7 @@ end
 use a +/- string to specify the occupation. If occbs=="-", the occupation from occas is used (closed-shell).
 if both are "-", the occupation is deduced from nelec.
 """
-function get_occvirt(EC::ECInfo, occas::String, occbs::String, norb, nelec)
+function get_occvirt(EC::ECInfo, occas::String, occbs::String, norb, nelec, ms2=0)
   if occas != "-"
     occa = parse_orbstring(occas)
     if occbs == "-"
@@ -83,8 +83,8 @@ function get_occvirt(EC::ECInfo, occas::String, occbs::String, norb, nelec)
       error("Inconsistency in OCCA ($occas) and OCCB ($occbs) definitions and the number of electrons ($nelec). Use ignore_error (-f) to ignore.")
     end
   else 
-    occa = [1:nelec÷2;]
-    occb = [1:(nelec+1)÷2;]
+    occa = [1:(nelec+ms2)÷2;]
+    occb = [1:(nelec-ms2)÷2;]
   end
   virta = [ i for i in 1:norb if i ∉ occa ]
   virtb = [ i for i in 1:norb if i ∉ occb ]
@@ -98,4 +98,4 @@ function get_occvirt(EC::ECInfo, occas::String, occbs::String, norb, nelec)
 end
 
 
-end #module
+end #module