Lots of molecules are working. 2s and 2p orbitals are running.

Author: willwheelera

GenerateStartingFunctions.py

@@ -28,7 +28,19 @@
 
 def Y1all(rvec, Z=1.0):
   rmag = np.sqrt(np.sum(rvec*rvec,1))
-  return np.sqrt(0.75/np.pi) * rvec/rmag
+  return np.sqrt(0.75/np.pi) * rvec / np.array([rmag,rmag,rmag]).T
+
+def Y1x(rvec, Z=1.0):
+  rmag = np.sqrt(np.sum(rvec*rvec,1))
+  return np.sqrt(0.75/np.pi) * rvec[:,0] / rmag
+
+def Y1y(rvec, Z=1.0):
+  rmag = np.sqrt(np.sum(rvec*rvec,1))
+  return np.sqrt(0.75/np.pi) * rvec[:,1] / rmag
+
+def Y1z(rvec, Z=1.0):
+  rmag = np.sqrt(np.sum(rvec*rvec,1))
+  return np.sqrt(0.75/np.pi) * rvec[:,2] / rmag
 
 # Radial functions
 # Note: this function requires that position coordinates are along dimension 1 of the array (not 0)
@@ -80,13 +92,13 @@ def psi_2p_all(e_pos_vec, i_pos, Z=1.0):
   return Y1all(e_pos_vec - i_pos) * R21(e_pos_vec - i_pos, Z)
 
 def psi_2px(e_pos_vec, i_pos, Z=1.0):
-  return Y1all(e_pos_vec - i_pos)[0] * R21(e_pos_vec - i_pos, Z)
+  return Y1x(e_pos_vec - i_pos)[0] * R21(e_pos_vec - i_pos, Z)
 
 def psi_2py(e_pos_vec, i_pos, Z=1.0):
-  return Y1all(e_pos_vec - i_pos)[1] * R21(e_pos_vec - i_pos, Z)
+  return Y1y(e_pos_vec - i_pos)[1] * R21(e_pos_vec - i_pos, Z)
 
 def psi_2pz(e_pos_vec, i_pos, Z=1.0):
-  return Y1all(e_pos_vec - i_pos)[2] * R21(e_pos_vec - i_pos, Z)
+  return Y1z(e_pos_vec - i_pos)[2] * R21(e_pos_vec - i_pos, Z)
 
 # Laplacian of S orbitals
 def Lpsi_1s(e_pos_vec, i_pos):
@@ -123,13 +135,16 @@ def psi_1s(self, e_pos_vec):
 class Atom:
   Z = 1.0
   i_pos = np.zeros(3)
-
+  
+  last_e_vec = np.zeros(3)
+  last_2p_vec = np.zeros(3)
+  
   def __init__(self, pos=np.array([0,0,0]), Z=1.0):
     self.i_pos = pos
     self.Z = float(Z)
 
-  last_e_vec = np.zeros(3)
-  last_2p_vec = np.zeros(3)
+  
+  
 
   def setPosition(self, pos):
     self.i_pos = pos
@@ -144,26 +159,17 @@ def psi_2p_all(self, e_pos_vec):
     return Y1all(e_pos_vec - self.i_pos) * R21(e_pos_vec - self.i_pos, self.Z)
 
   def psi_2px(self, e_pos_vec):
-    if (last_e_vec == e_pos_vec).all():
-      return last_2p_vec[0]
-    else:
-      last_e_vec = e_pos_vec.copy()
-      last_2p_vec = Y1all(e_pos_vec - self.i_pos)[0] * R21(e_pos_vec - self.i_pos, self.Z)
-      return last_2p_vec[0]
+    #if (self.last_e_vec == e_pos_vec).all():
+    #  return self.last_2p_vec[0]
+    #else:
+    #  self.last_e_vec = e_pos_vec.copy()
+    #  self.last_2p_vec = Y1all(e_pos_vec - self.i_pos) * R21(e_pos_vec - self.i_pos, self.Z)
+      return Y1x(e_pos_vec - self.i_pos) * R21(e_pos_vec - self.i_pos,self.Z)
 
   def psi_2py(self, e_pos_vec):
-    if (last_e_vec == e_pos_vec).all():
-      return last_2p_vec[1]
-    else:
-      last_e_vec = e_pos_vec.copy() 
-      last_2p_vec = Y1all(e_pos_vec - self.i_pos)[0] * R21(e_pos_vec - self.i_pos, self.Z)
-      return last_2p_vec[1]
+      return Y1y(e_pos_vec - self.i_pos) * R21(e_pos_vec - self.i_pos,self.Z)
+
 
   def psi_2pz(self, e_pos_vec):
-    if (last_e_vec == e_pos_vec).all():
-      return last_2p_vec[2]
-    else:
-      last_e_vec = e_pos_vec.copy() 
-      last_2p_vec = Y1all(e_pos_vec - self.i_pos)[0] * R21(e_pos_vec - self.i_pos, self.Z)
-      return last_2p_vec[2]
+      return Y1z(e_pos_vec - self.i_pos) * R21(e_pos_vec - self.i_pos,self.Z)
 

TrialFunctions2.py

@@ -101,6 +101,34 @@ def H2Molecule(ion_sep):
     #print 'Simulating H2Molecule'
     return wf
 
+def H3Molecule(ion_sep):
+    # ion_sep is in atomic units of Bohr radius 
+    ion_positions = np.array([
+        [-0.5*ion_sep, 0, 0], 
+        [0.5*ion_sep, 0, 0], 
+        [0,0.5*ion_sep, 0]]) * GSF.a_B
+    H_atom1 = GSF.H_atom(pos=np.array(ion_positions[0]))#,Z=1.0)
+    H_atom2 = GSF.H_atom(pos=np.array(ion_positions[1]))#,Z=1.0)
+    H_atom3 = GSF.H_atom(pos=np.array(ion_positions[2]))#,Z=1.0)
+    psi_laplacian = []
+    # two options for 2 electrons --> 2(up and down):0 or 1:1  (up: down or up:up)
+    # using 1:1 and up for both for now  
+    psi_array_up = np.array([H_atom1.psi_1s,H_atom2.psi_1s])
+    psi_array_down = np.array([H_atom3.psi_1s])
+
+    wf = WaveFunctionClass()
+    wf.setUpWavefunctions(psi_array_up)
+    wf.setDownWavefunctions(psi_array_down)
+    wf.setAtomicLaplacians(psi_laplacian)
+    wf.setAtomList([H_atom1,H_atom2,H_atom3])
+    #wf.setIonPositions(ion_positions)
+    #wf.setIonCharges(ion_charges)
+    wf.setNumUp(len(psi_array_up))
+    wf.setNumDown(len(psi_array_down))
+    
+    #print 'Simulating H2Molecule'
+    return wf
+
 def He2Molecule(ion_sep):
     # ion_sep is in atomic units of Bohr radius 
     ion_positions = np.array([

VMC.py

@@ -68,8 +68,8 @@ def MC_loop(WF, steps=1000, sig=0.5):
             index += 1
 
         E[t] = WF.LocalEnergy(Psi)
-        printtime = 1000
-        if (t+1)%printtime == 0: print 'Finished step '+str(t+1)+str(WF.e_positions)
+        printtime = 250
+        if (t+1)%printtime == 0: print 'Finished step '+str(t+1) #str(WF.e_positions)
 
     print 'Final prob ratio',prob_ratio
     print('Acceptance Rate:',(moves_accepted/(N*steps)))
@@ -125,7 +125,7 @@ def parseArgs(args,x):
     steps_input, bond_distance, sigma, jastrowB, jastrowD = parseArgs(sys.argv,x)
 
     #WF = GTF.H2Molecule(bond_distance, N_e=1)
-    WF = GTF.H2Molecule(bond_distance)
+    WF = GTF.HFMolecule(bond_distance)
     #WF = GTF.LithiumAtom()
     #WF = GTF.HeliumAtom()
     WF.Bee_same = jastrowB