Skip to content

Updated scripts, small changes to applied_efp.rst #23

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Jacklaw441 merged 1 commit into from
Jan 27, 2025
Merged

Updated scripts, small changes to applied_efp.rst #23

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 6 additions & 3 deletions doc/applied_efp.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,14 @@ FMO: Applied flexible EFP
Overview
--------

As is the case with many photoactive proteins,computational methods struggle to reproduce
experimental spectra for the Fenna-Matthews-Olson complex (FMO). Work by
Computational methods struggle to reproduce experimental spectra in photoactive proteins.
One such protein is the the Fenna-Matthews-Olson complex (FMO). Work by
`Kim et al <https://pubs.acs.org/doi/full/10.1021/acs.jpclett.9b03486>`_ shows that
flexible QM/EFP can be applied to FMO to correctly generate computational results in
quantitative agreement to experimental spectra.
quantitative agreement to experimental spectra. In short, the solvatochromic environments
surrounding each bacteriochloropyll a (BChl a) pigment are heavily polarizable; EFP captures
this polarizability explicitly where a standard QM/MM calculation will onyl approximate it,
thus the polarizable environment is crucial to accurately model photosynthesis.

The key to applying EFP to your system is to carefully define the active site and EFP region.
FMO is a trimeric protein with eight bacteriochloropyll a (BChl) pigments in each monomer.
Expand Down
15 changes: 8 additions & 7 deletions examples/flex-EFP/Scripts/cut_qm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -243,20 +243,21 @@ def get_screen(lines,coords,title):
if title in line:
start=1

def get_header(lines):
def get_header(lines,input_name):
header=[]
fragname=input_name.split('.')[0]
for line in lines:
header.append(line)
header.append(line.replace('$FRAGNAME',fragname))
#print(line)
if 'COORDINATES (BOHR)' in line:
return header

def main(inp, efp):
with open(inp,'r') as inp:
inp_lines=inp.readlines()
with open(efp,'r') as efp:
efp_lines=efp.readlines()
header=get_header(efp_lines)
with open(inp,'r') as inp_:
inp_lines=inp_.readlines()
with open(efp,'r') as efp_:
efp_lines=efp_.readlines()
header=get_header(efp_lines,inp)
rem_atoms, rem_pols = get_specials(inp_lines)
keep_coords,rem_coords=get_coords(efp_lines,rem_atoms,rem_pols)
keep_monop=get_monopoles(efp_lines,keep_coords)
Expand Down
230 changes: 23 additions & 207 deletions examples/flex-EFP/Scripts/fragment_RMSD.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,74 +12,14 @@
ang_cutoff=0.20 # Minimum RMSD allowed for "good" match
cutoff=ang_cutoff*1.8897259886 # Cutoff convert to Bohr

inp=sys.argv[1] # a_22_304.inp for example. Input file with fragment coords
inp=sys.argv[1]
#inp='g_952.inp' # a_22_304.inp for example. Input file with fragment coords
with open(inp, "r") as orig:
orgs = orig.readlines()

outname=inp.replace('.inp','.efp') # Output file name, same as input but changed extension.
# If no match is found, then outname will not be used.

def rotate_vector(vector, rotation_matrix):
#Different from coords because dipoles assumed centered on coordinate. NOT placed in Cartesian space
return np.dot(vector,rotation_matrix)

def rotate_quadrupole(quadrupole, rotation_matrix):
#R x Q x R_t
#Quadrupoles are 3x3s. But .efp give only six values, order is assumed xx, xy, xz, yy, yz, zz
#With symmetric off-diagonal elements.
q = np.array([
[quadrupole[0], quadrupole[3], quadrupole[4]],
[quadrupole[3], quadrupole[1], quadrupole[5]],
[quadrupole[4], quadrupole[5], quadrupole[2]],
])
q_rotated = rotation_matrix.T @ q @ rotation_matrix
return [q_rotated[0, 0], q_rotated[1, 1], q_rotated[2, 2],
q_rotated[0, 1], q_rotated[0, 2], q_rotated[1, 2]]

def rotate_octupole(octupole, rotation_matrix):
#Oxyz=Rxi*Ryj*Rzk*Oijk summed over i=x,y,z;j=x,y,z;z=x,y,z;
#Each term Oxyz of new tensor is composed of
#27 pieces generated from old tensor and rotation matrix!
#The indexes are hard-coded for some better readiability
o = np.array([
[[octupole[0], octupole[3], octupole[4]], [octupole[3], octupole[5], octupole[9]], [octupole[4], octupole[9], octupole[7]]],
[[octupole[3], octupole[5], octupole[9]], [octupole[5], octupole[1], octupole[6]], [octupole[9], octupole[6], octupole[8]]],
[[octupole[4], octupole[9], octupole[7]], [octupole[9], octupole[6], octupole[8]], [octupole[7], octupole[8], octupole[2]]],
])
#fancy function below does what we want. Look up np.einsum if you need clarification.
#o_rotated = np.einsum('ij,jkl,lm->ikm', rotation_matrix, o, rotation_matrix.T)
o_rotated = np.zeros((3, 3, 3))
for p in range(0,3):
for q in range(0,3):
for r in range(0,3):
o_rotated[p][q][r]=single_term_oct(p,q,r,o,rotation_matrix.T)
return [
o_rotated[0, 0, 0], o_rotated[1, 1, 1], o_rotated[2, 2, 2],
o_rotated[0, 0, 1], o_rotated[0, 0, 2], o_rotated[0, 1, 1],
o_rotated[1, 1, 2], o_rotated[0, 2, 2], o_rotated[1, 2, 2],
o_rotated[0, 1, 2]
]

def single_term_oct(p,q,r,octupole,rot_mat):
term=0.0
for i in range(0,3):
for j in range(0,3):
for k in range(0,3):
term+=rot_mat[p][i]*rot_mat[q][j]*rot_mat[r][k]*octupole[i][j][k]
return term

def rotate_polarizability(polarizability, rotation_matrix):
#RxQxR_t
q = np.array([
[polarizability[0], polarizability[3], polarizability[4]],
[polarizability[6], polarizability[1], polarizability[5]],
[polarizability[7], polarizability[8], polarizability[2]]])
q_rotated = rotation_matrix.T @ q @ rotation_matrix
return [q_rotated[0, 0],q_rotated[1, 1],q_rotated[2, 2],
q_rotated[0, 1],q_rotated[0, 2],q_rotated[1, 2],
q_rotated[1, 0],q_rotated[2, 0],q_rotated[2, 1]]
#return q_rotated.flatten()

def get_RMSD(coords1,coords2,mass):
tot=0.0
length=len(coords1)
Expand All @@ -106,8 +46,8 @@ def kabsch_algorithm(coords, coords2):
u, _, vh = np.linalg.svd(covariance_matrix)
rotation_matrix = np.dot(u, vh)

#if np.linalg.det(rotation_matrix) < 0:
# rotation_matrix[:, -1] *= -1
if np.linalg.det(rotation_matrix) < 0:
rotation_matrix[:, -1] *= -1
return rotation_matrix, com1, com2

def apply_transform(coords, rotation_matrix, com1, com2):
Expand All @@ -118,28 +58,34 @@ def apply_transform(coords, rotation_matrix, com1, com2):
return rotated_coords

amino_acid_dict = {'a':'ala','r':'arg','n':'asn','d':'asp','c':'cys',
'q':'gln','e':'glu','g':'gly','h':'his','i':'ile',
'q':'gln','e':'glu','g':'gly','h':'hip','i':'ile',
'l':'leu','k':'lys','m':'met','f':'phe','p':'pro',
's':'ser','t':'thr','w':'trp','y':'tyr','v':'val',
'hp':'hip','hd':'hid','he':'hie'}
atom_weights = {'H':1.0,'O':15.9949100,'C':12.0000000,'N':12.0030700,'S':32.0650000,'M':23.3040000}

try:
res=amino_acid_dict[inp[0]]
res=amino_acid_dict[inp.split('_')[0]]
except:
print('No library folder for: '+inp)
#directory = '/depot/lslipche/data/yb_boss/flexible_efp/efpdb/'+res
exit()

directory+=res

file_extension = '.efp' # change this to your desired file extension


#Grab coordinates from the input file, convert from Angstroms to bohrs.
start=0
xyz=[]
orig_coords=[]
weights=[]
inp_at_lines=[]
for line in orgs:
if '$end' in line:
start=0
elif(start==1):
inp_at_lines.append(line)
xyz.append(float(line.split()[2])/0.529177249)
xyz.append(float(line.split()[3])/0.529177249)
xyz.append(float(line.split()[4])/0.529177249)
Expand All @@ -148,27 +94,23 @@ def apply_transform(coords, rotation_matrix, com1, com2):
weights.append(atom_weights[line.split()[0][0]])
elif 'C1' in line:
start=1

minrmsd=20.0
minrmsd=20.0 # Arbitrary starting value. This will be replaced with real RMSD values in the loop.

# Loop through files in the directory
for filename in os.listdir(directory):
if filename.endswith(file_extension):
full_path = os.path.join(directory, filename)
#full_path=directory+'/gly5817.efp'
#print(f"Processing file: {full_path}")
with open(full_path, "r") as term:
terms = term.readlines()

else:
continue
start=0
xyz=[]
term_coords=[]
for line in terms:
if 'BO21' in line:
start=0
elif(start==1):
#print(line.split()[1])
#print(float(line.split()[1])*0.529177249)
xyz.append(float(line.split()[1]))
xyz.append(float(line.split()[2]))
xyz.append(float(line.split()[3]))
Expand All @@ -181,139 +123,13 @@ def apply_transform(coords, rotation_matrix, com1, com2):
rotation_matrix, com1, com2 = kabsch_algorithm(new_coords, coords2)
aligned_new_coords = apply_transform(new_coords, rotation_matrix, com1,com2)
rmsd=(get_RMSD(aligned_new_coords,coords2,weights))
#print(rmsd)
if(rmsd<minrmsd):
minrmsd=rmsd
if(rmsd<cutoff):
match=full_path
#break
#print(minrmsd)
#%%

#If good match found, transform that .efp to coords in .inp
# Coords, dipoles, quadrupoles, octupoles, polarizable points and tensor
# Monopoles and screen/screen2 do not need changed

#print(rmsd/1.8897259886, full_path, inp)
#print(get_RMSD(aligned_new_coords,coords2,weights))
#print(aligned_new_coords)
#print(coords2)
if(rmsd<minrmsd):
minrmsd=rmsd
match=full_path

if(rmsd<cutoff):
with open(match, "r") as efp:
f0 = efp.readlines()
start=0
#coords=[]
temp=[]
for line in f0:
newline=line
if ' STOP' in line:
start=0
elif(start==1):
#coords.append(apply_transform([float(line.split()[1]),float(line.split()[2]),float(line.split()[3])],rotation_matrix,com1,com2))
coords=(apply_transform([float(line.split()[1]),float(line.split()[2]),float(line.split()[3])],rotation_matrix,com1,com2))
newline=line.split()[0]
while(len(newline)<7):
newline=newline+' '
for i in range(len(coords)):
newline+=str(f"{coords[i]:16.10f}")
#newline=newline.replace(line.split()[i+1],str(f"{coords[i]:16.10f}"))
newline+=str(f"{float(line.split()[4]):12.7f}")
newline+=str(f"{float(line.split()[5]):5.1f}")+'\n'
elif(start==2):
dipole=rotate_vector([float(line.split()[1]),float(line.split()[2]),float(line.split()[3])],rotation_matrix)
newline=line.split()[0]
while(len(newline)<7):
newline=newline+' '
for i in range(len(dipole)):
newline+=f"{dipole[i]:16.10f}"
newline+='\n'
elif(start==3):
newline=''
if '>' in line:
name=line.split()[0]
while(len(name)<7):
name=name+' '
quadrupole=[float(line.split()[1]),float(line.split()[2]),float(line.split()[3]),float(line.split()[4])]
else:
quadrupole.append(float(line.split()[0]))
quadrupole.append(float(line.split()[1]))
rotated_quad=rotate_quadrupole(quadrupole, rotation_matrix)
temp.append(name+' '+f"{rotated_quad[0]:14.10f}"+' '+f"{rotated_quad[1]:14.10f}"+' '
+f"{rotated_quad[2]:14.10f}"+' '+f"{rotated_quad[3]:14.10f}"+' >\n'+
' '+f"{rotated_quad[4]:14.10f}"+' '+f"{rotated_quad[5]:14.10f}"+'\n')
elif(start==4):
newline=''
if(j%3==0):
name=line.split()[0]
while(len(name)<7):
name=name+' '
#print(line)
octupole=[float(line.split()[1]),float(line.split()[2]),float(line.split()[3]),float(line.split()[4])]
j+=1
elif(j%3==1):
for i in range(4):
octupole.append(float(line.split()[i]))
j+=1
else:
octupole.append(float(line.split()[0]))
octupole.append(float(line.split()[1]))
rotated_oct=rotate_octupole(octupole, rotation_matrix)
#print(rotated_oct)
temp.append(name+' '+f"{rotated_oct[0]:13.9f}"+' '+f"{rotated_oct[1]:13.9f}"+' '
+f"{rotated_oct[2]:13.9f}"+' '+f"{rotated_oct[3]:13.9f}"+' >\n'+
' '+f"{rotated_oct[4]:13.9f}"+' '+f"{rotated_oct[5]:13.9f}"' '
+f"{rotated_oct[6]:13.9f}"+' '+f"{rotated_oct[7]:13.9f}"+' >\n'+
' '+f"{rotated_oct[8]:13.9f}"+' '+f"{rotated_oct[9]:13.9f}"+'\n')
j=0
elif(start==5):
if 'CT' in line:
newline=line.split()[0]
while(len(newline)<4):
newline=newline+' '
coords=(apply_transform([float(line.split()[1]),float(line.split()[2]),float(line.split()[3])],rotation_matrix,com1,com2))
for i in range(len(coords)):
newline+=f"{coords[i]:16.10f}"
newline+='\n'
temp.append(newline)
elif(j%3==0):
polar=[float(line.split()[0]),float(line.split()[1]),float(line.split()[2]),float(line.split()[3])]
j+=1
elif(j%3==1):
for i in range(4):
polar.append(float(line.split()[i]))
j+=1
else:
polar.append(float(line.split()[0]))
rotated_polar=rotate_polarizability(polar,rotation_matrix)
temp.append(' '+f"{rotated_polar[0]:14.10f}"+' '+f"{rotated_polar[1]:14.10f}"+' '+
f"{rotated_polar[2]:14.10f}"+' '+f"{rotated_polar[3]:14.10f}"+' >\n'+
' '+f"{rotated_polar[4]:14.10f}"+' '+f"{rotated_polar[5]:14.10f}"+' '+
f"{rotated_polar[6]:14.10f}"+' '+f"{rotated_polar[7]:14.10f}"+' >\n'+
f"{rotated_polar[8]:16.10f}"+'\n')
j=0
elif 'COORDINATES (BOHR)' in line:
start=1
elif ' DIPOLES' in line:
start=2
elif ' QUADRUPOLES' in line:
start=3
temp.append(line)
elif ' OCTUPOLES' in line:
temp.append(line)
start=4
j=0
elif ' POLARIZABLE POINTS' in line:
start=5
temp.append(line)
j=0
if(start<3):
temp.append(newline)

outfile=inp.split('.')[0]+'.efp'
with open(outfile,'w') as outfile:
for line1 in temp:
outfile.write(line1)
if(minrmsd<cutoff):
os.system("python step4.Flexible_V5.py "+full_path+' '+inp+' -d -xr')
else:
with open('no_matches.txt','a') as outfile2:
outfile2.write(inp+'\n')
print('no match, run GAMESS for: '+inp)
#os.system('gms_slurm -p 20 -q standby -v 2023 '+inp)
Loading
Loading