More robust handling of 2-coordinate species

Instead of using the rotational method, first do a sum of Euclidean vectors. If sum is > sum_tol (0.5 A), then set adsorbate at this position. If not, use rotational approach.

mai/ads_sites.py

@@ -182,13 +182,17 @@ def get_nonplanar_ads_pos(self,scaled_sum_dist,center_coord):
 
 		return ads_pos
 
-	def get_bi_ads_pos(self,normal_vec,center_coord,site_idx):
+	def get_bi_ads_pos(self,normal_vec,scaled_sum_dist,center_coord,site_idx):
 		"""
 		Get adsorption site for a 2-coordinate site
 
 		Args:
 			normal_vec (numpy array): 1D numpy array for the
 			normal vector to the line
+
+			scaled_sum_dist (numpy array): 2D numpy array for the
+			scaled Euclidean distance vectors between each coordinating
+			atom and the central atom (i.e. the adsorption site)
 			
 			center_coord (numpy array): 1D numpy array for adsorption site
 			
@@ -199,52 +203,61 @@ def get_bi_ads_pos(self,normal_vec,center_coord,site_idx):
 			adsorption position
 		"""
 		r_cut = self.r_cut
+		sum_tol = self.sum_tol
 		overlap_tol = self.overlap_tol
 		dist = self.get_dist_planar(normal_vec)
+		n_overlap = None
 
-		#Prepare 2 temporary adsorbates
-		try_angles = np.arange(0,360,10)
-		ads_pos_temp_unrotated1 = center_coord + dist
-		ads_pos_temp_unrotated2 = center_coord - dist
-		ads_temp1 = Atoms([Atom('X',ads_pos_temp_unrotated1)])
-		ads_temp2 = Atoms([Atom('X',ads_pos_temp_unrotated2)])
-		mof_temp_orig = self.start_atoms.copy()
-
-		#Rotate one of the adsorbates about the axis
-		#defined by the two coordinating atoms
-		for i, angle in enumerate(try_angles):
-
-			#Add 2 overlapping atoms
-			mof_temp = mof_temp_orig.copy()
-			mof_temp.extend(ads_temp1)
-			mof_temp.extend(ads_temp2)
-
-			#Offset one of the atoms by small value
-			#(just to be higher than machine epsilon)
-			mof_temp[-1].position += 1e-6
-
-			#Set the new angle
-			mof_temp.set_angle(-1,site_idx,-2,angle)
-
-			#Determine number of nearby atoms
-			dist_mat = mof_temp.get_distances(-2,np.arange(0,
-				len(mof_temp)-2).tolist(),mic=True)
-			NNs = sum(dist_mat <= r_cut)
-			n_overlap = sum(dist_mat <= overlap_tol)
-			min_dist = np.min(dist_mat)
-
-			#Select best option
-			if i == 0:
-				ads_pos = mof_temp[-2].position
-				old_min_NNs = NNs
-				old_min_dist = min_dist
-				old_overlap = n_overlap
-			elif n_overlap <= old_overlap and (NNs < old_min_NNs or
-				(NNs == old_min_NNs and min_dist > old_min_dist)):
-				ads_pos = mof_temp[-2].position
-				old_min_NNs = NNs
-				old_min_dist = min_dist
-				old_overlap = n_overlap
+		if np.linalg.norm(scaled_sum_dist) > sum_tol:
+			ads_pos_nonplanar = self.get_nonplanar_ads_pos(scaled_sum_dist,center_coord)
+			NN_nonplanar, min_dist_nonplanar, n_overlap = self.get_NNs(ads_pos_nonplanar,
+				site_idx)
+		if n_overlap == 0:
+			ads_pos = ads_pos_nonplanar
+		else:
+			#Prepare 2 temporary adsorbates
+			try_angles = np.arange(0,360,10)
+			ads_pos_temp_unrotated1 = center_coord + dist
+			ads_pos_temp_unrotated2 = center_coord - dist
+			ads_temp1 = Atoms([Atom('X',ads_pos_temp_unrotated1)])
+			ads_temp2 = Atoms([Atom('X',ads_pos_temp_unrotated2)])
+			mof_temp_orig = self.start_atoms.copy()
+
+			#Rotate one of the adsorbates about the axis
+			#defined by the two coordinating atoms
+			for i, angle in enumerate(try_angles):
+
+				#Add 2 overlapping atoms
+				mof_temp = mof_temp_orig.copy()
+				mof_temp.extend(ads_temp1)
+				mof_temp.extend(ads_temp2)
+
+				#Offset one of the atoms by small value
+				#(just to be higher than machine epsilon)
+				mof_temp[-1].position += 1e-6
+
+				#Set the new angle
+				mof_temp.set_angle(-1,site_idx,-2,angle)
+
+				#Determine number of nearby atoms
+				dist_mat = mof_temp.get_distances(-2,np.arange(0,
+					len(mof_temp)-2).tolist(),mic=True)
+				NNs = sum(dist_mat <= r_cut)
+				n_overlap = sum(dist_mat <= overlap_tol)
+				min_dist = np.min(dist_mat)
+
+				#Select best option
+				if i == 0:
+					ads_pos = mof_temp[-2].position
+					old_min_NNs = NNs
+					old_min_dist = min_dist
+					old_overlap = n_overlap
+				elif n_overlap <= old_overlap and (NNs < old_min_NNs or
+					(NNs == old_min_NNs and min_dist > old_min_dist)):
+					ads_pos = mof_temp[-2].position
+					old_min_NNs = NNs
+					old_min_dist = min_dist
+					old_overlap = n_overlap
 
 		return ads_pos
 
@@ -319,6 +332,12 @@ def get_opt_ads_pos(self,mic_coords,site_idx):
 		self.io_stats['cnum'] = str(cnum)
 
 		#Calculate relevant quantities based on coordination number
+		if cnum >= 2:
+			#Calculate sum of Euclidean vectors
+			scaled_mic_coords = mic_coords*scale_factor/np.linalg.norm(
+				mic_coords,axis=1)[np.newaxis].T
+			scaled_sum_dist = np.sum(scaled_mic_coords,axis=0)
+			norm_scaled = np.linalg.norm(scaled_sum_dist)
 		if cnum == 1:
 			#Define direction as colinear with center atom and
 			#coordinating atom
@@ -328,20 +347,16 @@ def get_opt_ads_pos(self,mic_coords,site_idx):
 			#the two coordinating atoms
 			normal_vec = OLS_fit(mic_coords)
 		elif cnum >= 3:
-			#Calculate normal vector to best-fit plane and the
-			#sum of Euclidean vectors
-			scaled_mic_coords = mic_coords*scale_factor/np.linalg.norm(
-				mic_coords,axis=1)[np.newaxis].T
-			scaled_sum_dist = np.sum(scaled_mic_coords,axis=0)
-			norm_scaled = np.linalg.norm(scaled_sum_dist)
+			#Calculate normal vector to best-fit plane
 			rmse, normal_vec = TLS_fit(mic_coords)
 
 		#Get adsorption site based on coordination number
 		if cnum == 1:
 			dist = self.get_dist_planar(normal_vec)
 			ads_pos = center_coord-dist
 		elif cnum == 2:
-			ads_pos = self.get_bi_ads_pos(normal_vec,center_coord,site_idx)
+			ads_pos = self.get_bi_ads_pos(normal_vec,scaled_sum_dist,center_coord,
+				site_idx)
 		elif cnum == 3 and norm_scaled > sum_tol:
 			ads_pos = self.get_tri_ads_pos(normal_vec,scaled_sum_dist,center_coord,
 				site_idx)

mai/adsorbate_constructor.py

@@ -275,10 +275,10 @@ def get_adsorbate_grid(self,atoms_path=None,site_idx=None,grid_path=None,
 
 		site_pos = atoms[site_idx].position
 		ads_pos = get_best_grid_pos(atoms,max_dist,site_idx,grid_filepath)
-		if ads_pos is 'nogrid':
+		if ads_pos == 'nogrid':
 			print('WARNING: no grid for '+name)
 			return None
-		elif ads_pos is 'invalid':
+		elif ads_pos == 'invalid':
 			print('WARNING: all NaNs within cutoff for '+name)
 			return None
 		ads_optimizer = ads_pos_optimizer(self,new_mofs_path=new_mofs_path,

setup.py

@@ -6,7 +6,7 @@
       author_email='[email protected]',
       url='https://github.com/arosen93/mof-adsorbate-initializer',
       requires_python='>=3.6.0',
-      version='1.1',
+      version='1.2',
       packages=find_packages(),
       license='MIT'
      )