#1. Go to your C drive, make a file called "Workbench" if it is not already made. #2. Place your output file in Workbench. # Make sure your file begins with Optim_ and ends with .z # #3. Go to line 135. You should see a small block of red text surrounded by # three quotation marks. Change this block of text to match the constraints # you will want for your next optimization step. # # In the last line of code: #4. Place the name of your file in the red quotation marks. # Keep the quotation marks. # Do not include the Optim_ or .z in "file_name" # Now fill out torsion_scan, angle_s, angle_r, torsion_fixed, theta_fixed, # n_isomers, charge and n_atoms in the last line of code. # Follow the format presented in the second-to-last line of code. If you are # wondering what any of these variables are, look at the definitions provided below. # Click the "Run" button (looks like a green triangle, in top left of screen # under the "Help" tab.) # # # DEFINITIONS: # # torsion_scan = The torsion you have already scanned across. # e.g. if you have scanned across psi and held phi fixed, psi is torsion_scan. # # angle_s = Starting angle for torsion_scan. # # angle_r = Step-size for torsion_scan. # e.g. if you scanned every 15 degrees along psi, angle_r is 15. # # torsion_fixed = The torsion you held fixed. # e.g. if you have scanned across psi and held phi fixed, phi is torsion_fixed. # # theta_fixed = The value at which you are fixing torsion_fixed. # e.g. if you fixed phi at 90 degrees, theta_fixed is 90. # # n_isomers = The number of conformers you expect to have in your output file. # e.g. if you are scanning across all psi with a step-size of 15, you # should have 24 conformers. e.g. 360 / 15 == 24 # # charge = The total charge on your molecule. # e.g. if you have a neutral molecule, you should have a charge of 0. # # n_atoms = The total number of atoms in your molecule. # # # Workbench should now have xyz files for all isomers and the input files you # will need for the next optimization step. # #5. Remove all files from Workbench. #6. Repeat for each file you want to process. def make_lst(query,lst,file_name,wd): Optim = open(str(wd)+'Optim_'+str(file_name)+'.z', 'r') k = 20 linenum = 0 num_query = 0 for line in Optim: line = line.rstrip() linenum = linenum + 1 for i in xrange(len(line)-k+1): if query == line[i : i+k]: num_query = num_query +1 lst[num_query] = lst.get(str(num_query), linenum) Optim.close() def raise_the_alarm(lst,torsion_scan,angle_s,angle_r,torsion_fixed,theta_fixed,file_name,wd): Optim = open(str(wd)+'Optim_'+str(file_name)+'.z', 'r') query = "Non-Optimized Parameters" alarm = {} k = 24 linenum = 0 num_query = 0 for line in Optim: line = line.rstrip() linenum = linenum + 1 for i in xrange(len(line)-k+1): if query == line[i : i+k]: num_query = num_query +1 alarm[num_query] = alarm.get(str(num_query), linenum) if alarm == {}: print "All isomers were successfully optimized!" else: for item in lst: if alarm[num_query] == lst[item]: theta = angle_r*int(item) - angle_r + angle_s if theta >= 360: theta_scan = theta - 360 elif theta < 0: theta_scan = theta + 360 else: theta_scan = theta print "WARNING: The isomer at "+str(torsion_scan)+str(theta_scan)+"_"+str(torsion_fixed)+str(theta_fixed)+" did not optimize!" Optim.close() def compare_lists(lst1,lst2,output): for i in xrange(len(lst1)): for stndrd in lst2: if lst2[stndrd] <= lst1[i]: output[i] = lst2[stndrd] def label_lines(lst,file_name,wd): Optim = open(str(wd)+'Optim_'+str(file_name)+'.z', 'r') linenum = 0 for line in Optim: line = line.rstrip() lst[linenum] = lst.get(linenum,str(line)) linenum = linenum + 1 Optim.close() def xyz_writer(lst,torsion_scan,theta_scan,torsion_fixed,theta_fixed,n_atoms,wd): xyz = open(str(wd)+'Coord_'+str(torsion_scan)+str(theta_scan)+'_'+str(torsion_fixed)+str(theta_fixed)+'.xyz','w') xyz.write(str(n_atoms)+'\n') xyz.write("Coord_"+str(torsion_scan)+str(theta_scan)+'_'+str(torsion_fixed)+str(theta_fixed)+'\n') for i in xrange(1,n_atoms+1): for item in lst: if lst[item] == i: xyz.write(str(item)+'\n') xyz.close() def optim_writer(lst,torsion_scan,theta_scan,torsion_fixed,theta_fixed,charge,n_atoms,direction,wd): optim = open(str(wd)+'Optim_'+str(torsion_scan)+str(theta_scan)+'_'+str(torsion_fixed)+str(theta_fixed)+'_'+str(direction)+'.i','w') optim.write('%nprocshared=12'+'\n') optim.write('%'+'lindaworker=localhost'+'\n') optim.write('%mem=2GB'+'\n') optim.write('%' +'chk=Optim_'+str(torsion_scan)+str(theta_scan)+'_'+str(torsion_fixed)+str(theta_fixed)+'_'+str(direction)+'.chk'+'\n') optim.write('# b3lyp/6-31G* opt=modredun int=grid=ultrafine scrf=iefpcm'+'\n'+'\n') optim.write('Scan '+str(torsion_fixed)+' starting from '+str(theta_fixed)+' with '+str(torsion_scan)+' fixed at '+str(theta_scan)+'\n'+'\n') optim.write(str(charge)+' 1'+'\n') for a in xrange(1,n_atoms+1): for item in lst: if lst[item] == a: optim.write(str(item)+'\n') optim.write('\n') optim.write('D 2 5 7 19 S 23') if direction == 'f': optim.write(' 15.000000') if direction == 'r': optim.write(' -15.000000') optim.write(""" D 1 2 5 7 F D 7 18 29 37 F D 18 29 37 43 F D 7 19 33 44 F D 19 32 34 35 F D 32 36 38 46 F D 2 4 14 26 F D 4 6 15 27 F D 6 9 22 25 F D 9 3 8 10 F D 3 8 10 24 F """) #This is where you specify which angles are being fixed and which are being scanned. #The torsion that was fixed previously should replace the line "optim.write('D 2 5 7 19 S 23')" above. #D indicates Dihedral, the four numbers (each once space apart) indicate the atom numbers from the molecule, S indicates scan, and 23 is the number of isomers it is creating #Intuition may lead one to think the it should be 24 rather than 23 if the scan angle is 15 degrees but one angle has already been accounted for in the first scan #So this number comes from 360/(angle of rotation) - 1. optim.write('\n'+'\n') optim.close() def atom_convert(lst,output): for item in lst: converter = {1:"H",6:"C",7:"N",8:"O"} atom = item[17] new_atom = converter[int(atom)] new_item = str(item)[30:] new_item = str(new_atom) + str(new_item) output[new_item] = lst.get(str(new_item), lst[item]) def extract_isomer(inpt,key,torsion_scan,angle_s,angle_r,torsion_fixed,theta_fixed,n_isomers,charge,n_atoms): for i in xrange(1,n_isomers+1): lst_int = {} for line in inpt: if line > key[i]+3 and line < key[i]+n_atoms+4: lst_int[inpt[line]] = lst_int.get(str(inpt[line]),line-key[i]-3) lst_int2 = {} atom_convert(lst_int,lst_int2) theta = angle_r*i - angle_r + angle_s if theta > 360: theta_scan = theta - 360 elif theta <= 0: theta_scan = theta + 360 else: theta_scan = theta xyz_writer(lst_int2,torsion_scan,theta_scan,torsion_fixed,theta_fixed,n_atoms,wd) direction = 'f' optim_writer(lst_int2,torsion_scan,theta_scan,torsion_fixed,theta_fixed,charge,n_atoms,direction,wd) direction = 'r' optim_writer(lst_int2,torsion_scan,theta_scan,torsion_fixed,theta_fixed,charge,n_atoms,direction,wd) def run_code(file_name,torsion_scan,angle_s,angle_r,torsion_fixed,theta_fixed,n_isomers,charge,n_atoms,wd): qry_optim = "Optimized Parameters" lst_optim = {0:0} qry_stndrd = "Standard orientation" lst_stndrd = {0:0} make_lst(qry_optim,lst_optim,file_name,wd) raise_the_alarm(lst_optim,torsion_scan,angle_s,angle_r,torsion_fixed,theta_fixed,file_name,wd) make_lst(qry_stndrd,lst_stndrd,file_name,wd) lst_isomers = {} lines = {} compare_lists(lst_optim,lst_stndrd,lst_isomers) label_lines(lines,file_name,wd) extract_isomer(lines,lst_isomers,torsion_scan,angle_s,angle_r,torsion_fixed,theta_fixed,n_isomers,charge,n_atoms) #This is where you set your working directory! wd = "/Users/Desktop/Workbench/" #wd = "C:\\Workbench\\" #The format in python for your working directory is different depending on whether you are working on a mac or PC #FORMAT: run_code(file_name,torsion_scan,angle_s,angle_r,torsion_fixed,theta_fixed,n_isomers,charge,n_atoms) run_code("phi105_psi135","phi",105,15,"psi",135,24,0,48,wd)