Advertisement
Guest User

Untitled

a guest
Apr 25th, 2017
113
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 7.53 KB | None | 0 0
  1. from noodles import (gather, schedule)
  2. from qmworks import (adf, dftb, molkit, run, templates, Settings)
  3.  
  4. import numpy as np
  5.  
  6.  
  7. def main():
  8. backbone = '{}C1=CC=C(C=C1)C#CC2=CC=C({})C=C2'
  9.  
  10. # smiles of the electronDonors
  11. # electronDonors = ['N', 'P', 'S', 'CO']
  12. electronDonors = ['N']
  13. # smiles of the electronAcceptors
  14. # electronAcceptors = ['B', 'C(O)=O', 'N(=O)=O', 'C(F)(F)(F)']
  15. electronAcceptors = ['O']
  16.  
  17. # Genereate all the molecules
  18. molecules, names = create_molecules(backbone, electronDonors, electronAcceptors)
  19.  
  20. # Optimize with DFTB
  21. dftb_jobs = [dftb(templates.geometry, m, job_name='dftb_{}'.format(n)) for m, n in zip(molecules, names)]
  22.  
  23. # Schedule functions
  24. schedule_if_succesful = schedule(run_if_succesful)
  25. schedule_if_failure = schedule(run_if_failure)
  26. schedule_choose = schedule(choose_job)
  27.  
  28. # Optimize with DFT in gas phase
  29. s = create_settings_dft_gas()
  30. gas_jobs = [schedule_if_failure(adf, s, job, mol, name, prefix='adf_')
  31. for job, mol, name in zip(dftb_jobs, molecules, names)]
  32.  
  33. preoptimized_jobs = [schedule_choose(j1, j2) for j1, j2 in zip(dftb_jobs, gas_jobs)]
  34.  
  35. #Optimize with COSMO and compute frequency
  36. s = create_settings_dft_cosmo_freq()
  37. adf_cosmo_freq_jobs = [schedule_if_succesful(adf, s, job, name, prefix='solvation')
  38. for job, name in zip(preoptimized_jobs, names)]
  39.  
  40. # Filter jobs which frequencies are greater than zero
  41. schedule_check_freq = schedule(check_frequencies)
  42. candidates = [schedule_check_freq(j) for j in adf_cosmo_freq_jobs]
  43.  
  44. # Schedule a TD-DFT job with non-negative frequencies
  45. s = create_settings_TD_DFT()
  46. td_dft_jobs = [adf(s, job, job_name='td_dft_'.format(job.job_name)) for job in candidates]
  47.  
  48. schedule_getattr = schedule(extract_property)
  49. energies = [schedule_getattr(j, 'energy') for j in td_dft_jobs]
  50.  
  51. rs = run(gather(*energies), folder='biphenyl')
  52.  
  53. print(rs)
  54.  
  55. def create_molecules(backbone, electronDonors, electronAcceptors):
  56. """
  57. Create all the possible combination by changing one para substituent
  58. by an electronAcceptors and the other para position by an electronDonors
  59. """
  60. # Smiles representing the molecules to simulate
  61. smiles = [backbone.format(donor, acceptor)
  62. for donor in electronDonors for acceptor in electronAcceptors]
  63. # Plams molecules
  64. geometries = [molkit.from_smiles(s, forcefield='uff') for s in smiles]
  65.  
  66. # Name of the molecules
  67. invalid_chars = ['=', '(', ')', '[', ']', '#']
  68. names = [''.join(list(filter(lambda c: c not in invalid_chars, name))) for name in smiles]
  69.  
  70. return geometries, names
  71.  
  72.  
  73.  
  74. def adf_general_settings():
  75. """
  76. ADF settings use for the calculations
  77. """
  78. s = Settings()
  79.  
  80. # Basis
  81. s.specific.adf.basis.createoutput = None
  82.  
  83. # SCF
  84. s.specific.adf.SCF.iterations = 200
  85. s.specific.adf.SCF.convergence = '1.0e-6 1.0e-5'
  86.  
  87. # zlmfit
  88. s.specific.adf.zlmfit.quality = 'verygood'
  89.  
  90. # BeckeGrid
  91. s.specific.adf.BeckeGrid.quality = 'verygood'
  92.  
  93. # realtivistic
  94. s.specific.adf.Relativistic = 'SCALAR ZORA'
  95.  
  96. # Solvation
  97. with open('radii', 'r') as f:
  98. radii = f.read()
  99.  
  100. s.specific.adf.Solvation.div = 'ndiv=3 min=0.5 ofac=0.8'
  101. s.specific.adf.Solvation.solvent = 'name=Dichloromethane emp=0.0'
  102. s.specific.adf.Solvation.charged = 'conv=1e-10 iter=1000 corr'
  103. s.specific.adf.Solvation.CSMRSP = ''
  104. s.specific.adf.Solvation.Radii = '\n{}'.format(radii)
  105.  
  106. return s
  107.  
  108. def create_settings_TD_DFT(basis='TZ2P'):
  109. """
  110. create input for a TD-DFT calculation
  111. """
  112. s = adf_general_settings()
  113.  
  114. # Basis
  115. s.specific.adf.basis.type = '{}'.format(basis)
  116. s.specific.adf.basis.core = 'none'
  117.  
  118.  
  119. #COSMO
  120. s.specific.adf.Solvation.solvent = 'name=Dichloromethane emp=0.0 neql=2.028'
  121.  
  122. # functional
  123. s.specific.adf.xc.hybrid = 'CAMYB3LYP'
  124. s.specific.adf.xc.xcfun = ''
  125. s.specific.adf.xc.rangesep = ''
  126.  
  127. # excitations
  128. s.specific.adf.excitations.Davidson = ''
  129. s.specific.adf.excitations.lowest = 7
  130. s.specific.adf.eprint.sfo = "eig ovl"
  131.  
  132. return s
  133.  
  134. def create_settings_dft_gas(basis='DZP'):
  135. """
  136. Create input for an optimizatiion with ADF in gas phase
  137. """
  138. s = Settings()
  139.  
  140. # Basis
  141. s.specific.adf.basis.createoutput = None
  142.  
  143. # SCF
  144. s.specific.adf.SCF.iterations = 200
  145. s.specific.adf.SCF.convergence = '1.0e-6 1.0e-5'
  146.  
  147. # zlmfit
  148. s.specific.adf.zlmfit.quality = 'good'
  149.  
  150. # BeckeGrid
  151. s.specific.adf.BeckeGrid.quality = 'good'
  152.  
  153. # realtivistic
  154. s.specific.adf.Relativistic = 'SCALAR ZORA'
  155.  
  156.  
  157. # Basis
  158. s.specific.adf.basis.type = '{}'.format(basis)
  159. s.specific.adf.basis.core = 'small'
  160.  
  161. # functional
  162. s.specific.adf.xc.gga = 'BP86'
  163.  
  164. # geometry
  165. s.specific.adf.geometry.optim = 'delocalized'
  166. s.specific.adf.geometry.iterations = 99
  167. s.specific.adf.geometry.convergence = 'e=0.00001 grad=0.0001'
  168.  
  169. # SCF
  170. s.specific.adf.SCF.mix = 0.2
  171.  
  172. return s
  173.  
  174.  
  175.  
  176. def create_settings_dft_cosmo_freq(basis='TZ2P'):
  177. """
  178. Create input for an optimization and posterior frequencies calculation with ADF
  179. """
  180. s = adf_general_settings()
  181. s.specific.adf.analyticalfreq
  182. s.specific.adf.scanfreq = '-500 100'
  183. #COSMO
  184. s.specific.adf.Solvation.solvent = 'name=Dichloromethane emp=0.0'
  185.  
  186. # Basis
  187. s.specific.adf.basis.type = '{}'.format(basis)
  188. s.specific.adf.basis.core = 'small'
  189.  
  190. # functional
  191. s.specific.adf.xc.gga = 'BP86'
  192.  
  193. # geometry
  194. s.specific.adf.geometry.optim = 'delocalized'
  195. s.specific.adf.geometry.iterations = 99
  196. s.specific.adf.geometry.convergence = 'e=0.000001 grad=0.00001'
  197.  
  198. # SCF
  199. s.specific.adf.SCF.mix = 0.2
  200.  
  201. return s
  202.  
  203. def filter_None(jobs):
  204. """
  205. Remove the non-valid jobs
  206. """
  207. return list(filter(lambda j: j is not None, jobs))
  208.  
  209.  
  210. def check_frequencies(job):
  211. """
  212. Filter the jobs which frequencies are non-negative
  213. """
  214. predicate = np.any(np.array(job.frequencies) < 0)
  215. if predicate:
  216. print("Job: ", job.job_name, " has at least a negative frequency")
  217. return None
  218. else:
  219. return job
  220.  
  221.  
  222. def extract_property(job, prop):
  223. """
  224. wrapper around getattr
  225. """
  226. if job is None:
  227. return None
  228. else:
  229. return getattr(job, prop)
  230.  
  231.  
  232. def choose_job(job1, job2):
  233. """
  234. Try to get an attribute from job1 otherwise from job2 else None
  235. """
  236. if job1 is not None and job1.status not in ['failed', 'crashed']:
  237. return job1
  238. elif job2 is not None and job2.status not in ['failed', 'crashed']:
  239. return job2
  240. else:
  241. return None
  242.  
  243.  
  244. def run_if_succesful(fun, settings, job, name, prefix='job'):
  245. """
  246. Schedule a job if its depencies finished succesfully
  247. """
  248. if job is not None and job.status not in ['failed', 'crashed']:
  249. job_name = '{}_{}'.format(prefix, name)
  250. return fun(settings, job.molecule, job_name=job_name)
  251. else:
  252. return None
  253.  
  254. def run_if_failure(fun, settings, job, mol, name, prefix='job'):
  255. """
  256. Schedule a job if its depencies finished succesfully
  257. """
  258. if job is None or job.status in ['failed', 'crashed']:
  259. job_name = '{}_{}'.format(prefix, name)
  260. return fun(settings, mol, job_name=job_name)
  261. else:
  262. return None
  263.  
  264. if __name__ == "__main__":
  265. main()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement