Added back ZORA relativity; untested; removed adp_vector & ADP_tensor…

… sanity check
dylan-jayatilaka · Apr 23, 2024 · fa66718 · fa66718
1 parent 931739f
commit fa66718
Show file tree

Hide file tree

Showing 3 changed files with 129 additions and 126 deletions.
diff --git a/foofiles/crystal.foo b/foofiles/crystal.foo
@@ -5683,10 +5683,11 @@ contains
          ! Get ADP2s
          .fragment_atom(f).put_ADP2_vector_to(adp2)
 
+
          ! Sanity check
          ENSURE(NOT (refine_Uiso AND refine_ADP4s),"can't refine isotropic ADPs & ADP4s for "//trim(tag))
          ENSURE(NOT (refine_Uiso AND refine_ADP3s),"can't refine isotropic ADPs & ADP3s for "//trim(tag))
-         ENSURE(adp2.equals(.fragment_atom(f).ADP_tensor),"ADP2 not equal to ADP_tensor")
+       ! ENSURE(adp2.equals(.fragment_atom(f).ADP_tensor),"ADP2 not equal to ADP_tensor")
 
          ! ADP4's !!!!!!!!!!!!!
          if (refine_ADP4s) then
@@ -6148,7 +6149,7 @@ contains
          ! Sanity check
          ENSURE(NOT (refine_Uiso AND refine_ADP4s),"can't refine isotropic ADPs & ADP4s for "//trim(tag))
          ENSURE(NOT (refine_Uiso AND refine_ADP3s),"can't refine isotropic ADPs & ADP3s for "//trim(tag))
-         ENSURE(adp2.equals(.fragment_atom(f).ADP_tensor),"ADP2 not equal to ADP_tensor")
+       ! ENSURE(adp2.equals(.fragment_atom(f).ADP_tensor),"ADP2 not equal to ADP_tensor")
 
          ! ADP4's !!!!!!!!!!!!!
          if (refine_ADP4s) then
@@ -6516,7 +6517,7 @@ contains
       .fragment_atom(f).put_ADP2_vector_to(adp2)
 
       ! Sanity check
-      ENSURE(adp2.equals(.fragment_atom(f).ADP_tensor),"ADP2 not equal to ADP_tensor")
+    ! ENSURE(adp2.equals(.fragment_atom(f).ADP_tensor),"ADP2 not equal to ADP_tensor")
 
       ! Doing anharmomic terms?
       do_ADP4s = .fragment_atom(f).has_only_ADP4s_and_errors

diff --git a/foofiles/molecule.ints.foo b/foofiles/molecule.ints.foo
@@ -1280,93 +1280,93 @@ contains
 
    end
 
-!   make_1e_ZORA_matrices(T,Zx,Zy,Zz) ::: PURE
-!   ! Calculate the one-electron ZORA spin orbit matrices numerically.
-!   ! This includes the relativitically modified kinetic energy integrals.
-!      self :: INOUT
-!      T, Zx,Zy,Zz :: MAT{REAL}, OUT
-!
-!   ENSURE(.basis_info_made, "no basis info")
-!   ENSURE(.atom.allocated,  "no atom list")
-!   ENSURE(.becke_grid.allocated,"no becke_grid")
-!   ENSURE(.overlapping_atoms.allocated,"no overlapping_atoms")
-!
-!      ZORA  :: MAT4{REAL}@
-!      a1,b1 :: MAT3{REAL}@
-!      pt    :: MAT{REAL}@
-!      wt,V0 :: VEC{REAL}@
-!      atom1 :: VEC{INT}(1)
-!      atom2 :: VEC{INT}(2)
-!      atoms :: VEC{INT}@
-!      fa,la,na, a,aa,ma, atom_a :: INT
-!      fb,lb,nb, b,bb,mb, atom_b :: INT
-!      n_pt,q,k,l :: INT
-!      same_atoms :: BIN
-!
-!      ZORA.create(.n_bf,.n_bf,3,3)
-!
-!      ! Make the lower half of the ZORA spin orbit integrals
-!      do q = 1,.INQ:no_of_atom_pairs
-!
-!         if (NOT .overlapping_atoms(q)) cycle
-!
-!         .SET:set_atom_pair_indices_from(q,atom_a,atom_b,fa,la,na,fb,lb,nb)
-!
-!         same_atoms = atom_a==atom_b
-!         if (same_atoms) then; atom1 = [atom_a];        atoms = atom1
-!         else;                 atom2 = [atom_a,atom_b]; atoms = atom2
-!         end
-!
-!         ! Make the grid points and weights
-!         .becke_grid.make_grid(pt,wt,atoms,compress=TRUE,weight_is_0=FALSE)
-!         n_pt = pt.dim1                         ! Keep the weight_is_0 array
-!
-!         V0.create(n_pt)                        ! local potential grids
-!         .INTS:ZORA_potential(V0,pt)
-!         .becke_grid.make_bf_grids(a1,b1,ma,mb,pt,atoms)
-!         V0 = wt*V0
-!         wt.destroy; pt.destroy
-!         do k = 1,3
-!         do l = 1,3
-!         do a = 1,na
-!         do b = 1,nb
-!            aa = fa + a
-!            bb = fb + b
-!            ZORA(aa,bb,k,l) = sum(V0*a1(:,a,k)*b1(:,b,l))
-!         end
-!         end
-!         end
-!         end
-!         V0.destroy
-!
-!         b1.destroy
-!         a1.destroy
-!
-!      end
-!
-!      ! Make the upper half of the ZORA spin orbit integrals
-!      do k = 1,3
-!      do l = 1,3
-!      do a = 1,.n_bf
-!      do b = 1,a-1
-!         ZORA(b,a,l,k) = ZORA(a,b,k,l)
-!      end
-!      end
-!      end
-!      end
-!
-!      ! Assemble the ZORA contribution to the 1 electron hamiltonian
-!      ! Scalar kinetic energy contribution
-!      T  = ZORA(:,:,1,1) + ZORA(:,:,2,2) + ZORA(:,:,3,3)
-!
-!      ! Spin-dependent spin-orbit contribution
-!      Zx = ZORA(:,:,2,3) - ZORA(:,:,3,2)
-!      Zy = ZORA(:,:,3,1) - ZORA(:,:,1,3)
-!      Zz = ZORA(:,:,1,2) - ZORA(:,:,2,1)
-!
-!      ZORA.destroy
-!
-!   end
+   make_1e_ZORA_matrices(T,Zx,Zy,Zz) ::: PURE
+   ! Calculate the one-electron ZORA spin orbit matrices numerically.
+   ! This includes the relativitically modified kinetic energy integrals.
+      self :: INOUT
+      T, Zx,Zy,Zz :: MAT{REAL}, OUT
+
+   ENSURE(.basis_info_made, "no basis info")
+   ENSURE(.atom.allocated,  "no atom list")
+   ENSURE(.becke_grid.allocated,"no becke_grid")
+   ENSURE(.overlapping_atoms.allocated,"no overlapping_atoms")
+
+      ZORA  :: MAT4{REAL}@
+      a1,b1 :: MAT3{REAL}@
+      pt    :: MAT{REAL}@
+      wt,V0 :: VEC{REAL}@
+      atom1 :: VEC{INT}(1)
+      atom2 :: VEC{INT}(2)
+      atoms :: VEC{INT}@
+      fa,la,na, a,aa,ma, atom_a :: INT
+      fb,lb,nb, b,bb,mb, atom_b :: INT
+      n_pt,q,k,l :: INT
+      same_atoms :: BIN
+
+      ZORA.create(.n_bf,.n_bf,3,3)
+
+      ! Make the lower half of the ZORA spin orbit integrals
+      do q = 1,.INQ:no_of_atom_pairs
+
+         if (NOT .overlapping_atoms(q)) cycle
+
+         .SET:set_atom_pair_indices_from(q,atom_a,atom_b,fa,la,na,fb,lb,nb)
+
+         same_atoms = atom_a==atom_b
+         if (same_atoms) then; atom1 = [atom_a];        atoms = atom1
+         else;                 atom2 = [atom_a,atom_b]; atoms = atom2
+         end
+
+         ! Make the grid points and weights
+         .becke_grid.make_grid(pt,wt,atoms,compress=TRUE,weight_is_0=FALSE)
+         n_pt = pt.dim1                         ! Keep the weight_is_0 array
+
+         V0.create(n_pt)                        ! local potential grids
+         .INTS:ZORA_potential(V0,pt)
+         .becke_grid.make_bf_grids(a1,b1,ma,mb,pt,atoms)
+         V0 = wt*V0
+         wt.destroy; pt.destroy
+         do k = 1,3
+         do l = 1,3
+         do a = 1,na
+         do b = 1,nb
+            aa = fa + a
+            bb = fb + b
+            ZORA(aa,bb,k,l) = sum(V0*a1(:,a,k)*b1(:,b,l))
+         end
+         end
+         end
+         end
+         V0.destroy
+
+         b1.destroy
+         a1.destroy
+
+      end
+
+      ! Make the upper half of the ZORA spin orbit integrals
+      do k = 1,3
+      do l = 1,3
+      do a = 1,.n_bf
+      do b = 1,a-1
+         ZORA(b,a,l,k) = ZORA(a,b,k,l)
+      end
+      end
+      end
+      end
+
+      ! Assemble the ZORA contribution to the 1 electron hamiltonian
+      ! Scalar kinetic energy contribution
+      T  = ZORA(:,:,1,1) + ZORA(:,:,2,2) + ZORA(:,:,3,3)
+
+      ! Spin-dependent spin-orbit contribution
+      Zx = ZORA(:,:,2,3) - ZORA(:,:,3,2)
+      Zy = ZORA(:,:,3,1) - ZORA(:,:,1,3)
+      Zz = ZORA(:,:,1,2) - ZORA(:,:,2,1)
+
+      ZORA.destroy
+
+   end
 
 !   make_ENA_mx(Z)
 !   ! Calculate the one-electron electron nuclear attraction matrix numerically.

diff --git a/foofiles/molecule.scf.foo b/foofiles/molecule.scf.foo
@@ -2503,7 +2503,7 @@ contains
       select case (.scfdata.relativity_kind)
     ! case ("douglas-kroll-hess"); .:add_gc_DKH_core_mx
     ! case ("dkh");                .:add_gc_DKH_core_mx
-    ! case ("zora");               .:add_gc_ZORA_core_mx
+      case ("zora");               .:add_gc_ZORA_core_mx
     ! case ("iotc");               .:set_gc_IOTC_core_mx
       case ("pauli");              .:add_gc_Pauli_core_mx
       case ("none");               .:add_gc_core_mx
@@ -2549,41 +2549,43 @@ contains
 !
 !   end
 
-!   add_gc_ZORA_core_mx
-!   ! Add the core hamiltonain to a general complex "F"
-!   ! NOTE: only lower half is made
-!   ENSURE(.scfdata.allocated,"no scfdata")
-!   ENSURE(.scfdata.using_1e_zora_term," no Pauli?")
-!   ENSURE(.core_mx.is_allocated_with_genre("gc"),"no core matrix")
-!
-!      T,Lx,Ly,Lz :: MAT{REAL}@
-!      fac :: REAL
-!      I :: CPX
-!
-!      I = (ZERO,ONE)
-!
-!      T.create(.n_bf,.n_bf)
-!      Lx.create(.n_bf,.n_bf)
-!      Ly.create(.n_bf,.n_bf)
-!      Lz.create(.n_bf,.n_bf)
-!      .INTS:make_1e_ZORA_matrices(T,Lx,Ly,Lz)
-!
-!      .core_mx.gc.aa_block_plus(T)
-!      .core_mx.gc.bb_block_plus(T)
-!
-!      fac = G_FACTOR/TWO
-!      fac = fac * .scfdata.sl_1e_factor
-!      .core_mx.gc.ba_block_plus(Lx,-fac*I)
-!      .core_mx.gc.ba_block_plus(Ly,fac)
-!      .core_mx.gc.aa_block_plus(Lz,-fac*I)
-!      .core_mx.gc.bb_block_plus(Lz,fac*I)
-!      
-!      T.destroy
-!      Lz.destroy
-!      Ly.destroy
-!      Lx.destroy
-!
-!   end
+   add_gc_ZORA_core_mx ::: PURE
+   ! Add the core hamiltonain to a general complex "F"
+   ! NOTE: only lower half is made
+      self :: INOUT
+
+   ENSURE(.scfdata.allocated,"no scfdata")
+   ENSURE(.scfdata.using_1e_zora_term," no Pauli?")
+   ENSURE(.core_mx.is_allocated_with_genre("gc"),"no core matrix")
+
+      T,Lx,Ly,Lz :: MAT{REAL}@
+      fac :: REAL
+      I :: CPX
+
+      I = (ZERO,ONE)
+
+      T.create(.n_bf,.n_bf)
+      Lx.create(.n_bf,.n_bf)
+      Ly.create(.n_bf,.n_bf)
+      Lz.create(.n_bf,.n_bf)
+      .INTS:make_1e_ZORA_matrices(T,Lx,Ly,Lz)
+
+      .core_mx.gc.aa_block_plus(T)
+      .core_mx.gc.bb_block_plus(T)
+
+      fac = G_FACTOR/TWO
+      fac = fac * .scfdata.sl_1e_factor
+      .core_mx.gc.ba_block_plus(Lx,-fac*I)
+      .core_mx.gc.ba_block_plus(Ly, fac  )
+      .core_mx.gc.aa_block_plus(Lz,-fac*I)
+      .core_mx.gc.bb_block_plus(Lz, fac*I)
+
+      T.destroy
+      Lz.destroy
+      Ly.destroy
+      Lx.destroy
+
+   end
 
 !   set_gc_IOTC_core_mx
 !   ! Add the core hamiltonain to a general complex "F"