diff --git a/.github/workflows/docs_compile.yml b/.github/workflows/docs_compile.yml index 00e4583..664a292 100644 --- a/.github/workflows/docs_compile.yml +++ b/.github/workflows/docs_compile.yml @@ -30,28 +30,17 @@ jobs: large-packages: true docker-images: true swap-storage: true - - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v3 + - name: Setup micromamba + uses: mamba-org/setup-micromamba@v1 with: - python-version: ${{ matrix.python-version }} - - name: Install Miniconda - uses: conda-incubator/setup-miniconda@v2 - with: - auto-update-conda: true - python-version: ${{ matrix.python-version }} - use-mamba: true - - - name: Initialize conda - shell: bash -l {0} - run: conda init bash + environment-file: QCflow.yml + environment-name: QCflow + cache-environment: true - - name: Activate conda environment, install dependencies + - name: Install dependencies shell: bash -l {0} run: | - source ~/.bashrc - conda env create -f QCflow.yml - conda activate QCflow - conda install sphinx=8.2.3 + micromamba install sphinx=8.2.3 pip install sphinx_rtd_theme>=2.0.0 pip install standard-imghdr pip install . @@ -59,8 +48,6 @@ jobs: - name: Run build shell: bash -l {0} run: | - source ~/.bashrc - conda activate QCflow cd docs/ make html @@ -69,4 +56,4 @@ jobs: with: publish_branch: docs github_token: ${{ secrets.qcflow2}} - publish_dir: ./docs/_build/html \ No newline at end of file + publish_dir: ./docs/_build/html diff --git a/.github/workflows/run_test.yml b/.github/workflows/run_test.yml index 0759137..7aa1fa7 100644 --- a/.github/workflows/run_test.yml +++ b/.github/workflows/run_test.yml @@ -11,29 +11,23 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v2 - - name: Set up Conda - uses: conda-incubator/setup-miniconda@v2 + - name: Setup micromamba + uses: mamba-org/setup-micromamba@v1 with: - auto-activate-base: false - - name: Create and activate environment + environment-file: QCflow.yml + environment-name: QCflow + cache-environment: true + - name: Install dependencies + shell: bash -l {0} run: | - conda config --add channels defaults - conda config --add channels conda-forge - conda config --add channels anaconda - conda env create -f QCflow.yml - conda init bash - source ~/.bashrc - conda activate QCflow python -m pip install --upgrade pip - python -m pip install pytest # Ensure pytest is installed + python -m pip install pytest - name: Debug Environment + shell: bash -l {0} run: | - source ~/.bashrc - conda activate QCflow python --version python -m pip show pytest - name: Run Tests + shell: bash -l {0} run: | - source ~/.bashrc - conda activate QCflow python -m pytest -v ./tests/test.py diff --git a/Examples/FilteredSmilesFragments_props.csv b/Examples/FilteredSmilesFragments_props.csv new file mode 100644 index 0000000..51abd0e --- /dev/null +++ b/Examples/FilteredSmilesFragments_props.csv @@ -0,0 +1,111 @@ +smiles,valid,mw,n_rot,n_arom,psuedo_groups,fingerprint,subset +O=C(O)C=CC1=CC=C(O)C=C1,True,164.16,1,1,COC(=O)C=Cc1ccc(O)cc1,[0. 0. 0. ... 0. 0. 0.],mono +CC1=C(Cl)C=CC=C1NC1=NC=CC=C1C(=O)O,True,262.69599999999997,3,2,COC(=O)c1cccnc1Nc1cccc(C)c1C,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)C1=CC(I)=CC(I)=C1O,True,389.914,1,1,COC(=O)c1cc(C)cc(C)c1O,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)C=CC1=CC=CC=C1O,True,164.16,1,1,COC(=O)C=Cc1ccccc1O,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)C1=CNC2=CC=CC=C12,True,161.16,1,2,COC(=O)c1cn(C)c2ccccc12,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)CC1=CNC2=CC=CC=C12,True,175.187,2,2,COC(=O)Cc1cn(C)c2ccccc12,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)C=CC1=CC=CC(O)=C1,True,164.16,1,1,COC(=O)C=Cc1cccc(O)c1,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)[C@H]1O[C@@H](OC2=CC=C(C3=COC4=CC(O)=CC=C4C3=O)C=C2)[C@H](O)[C@@H](O)[C@@H]1O,True,430.3650000000001,3,3,COC(=O)[C@H]1O[C@@H](Oc2ccc(-c3coc4cc(O)ccc4c3=O)cc2)[C@H](O)[C@@H](O)[C@@H]1O,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)[C@H]1O[C@@H](OC2=CC=C3C(=O)C(C4=CC=C(O)C=C4)=COC3=C2)[C@H](O)[C@@H](O)[C@@H]1O,True,430.36500000000007,3,3,COC(=O)[C@H]1O[C@@H](Oc2ccc3c(=O)c(-c4ccc(O)cc4)coc3c2)[C@H](O)[C@@H](O)[C@@H]1O,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)C1=C(C2=CNC3=CC=CC=C23)C(C2=CNC3=CC=CC=C23)=C(C(=O)O)N1,True,385.37900000000013,4,5,COC(=O)c1c(-c2cn(C)c3ccccc23)c(-c2cn(C)c3ccccc23)c(C(=O)OC)n1C,[0. 0. 0. ... 0. 0. 0.],mono +O=C1NC(=O)C(C2=CNC3=CC=CC=C23)=C1C1=CNC2=CC=CC=C12,True,327.3430000000001,2,4,Cn1cc(C2=C(c3cn(C)c4ccccc34)C(=O)NC2=O)c2ccccc21,[0. 0. 0. ... 0. 0. 0.],mono +O=C1NCC2=C1C1=C(NC3=CC=CC=C31)C1=C2C2=CC=CC=C2N1,True,311.344,0,5,Cn1c2ccccc2c2c3c(c4c5ccccc5n(C)c4c21)C(=O)NC3,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)CC1=CNC2=CC=C(O)C=C12,True,191.186,2,2,COC(=O)Cc1cn(C)c2ccc(O)cc12,[0. 0. 0. ... 0. 0. 0.],mono +COC1=CC=C2C(=C1)C(CC(=O)O)=C(C)N2C(=O)C1=CC=C(Cl)C=C1,True,357.793,5,3,COC(=O)Cc1c(C)n(C(=O)c2ccc(C)cc2)c2ccc(OC)cc12,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)CCC1=NC=CC2=C1NC1=CC=CC=C12,True,240.26200000000003,3,3,COC(=O)CCc1nccc2c3ccccc3n(C)c12,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)C1=CC2=CC=CC=C2N1,True,161.16,1,2,COC(=O)c1cc2ccccc2n1C,[0. 0. 0. ... 0. 0. 0.],mono +COC1=CC=C2C(=C1)NC1=C(CCC(=O)O)N=CC=C12,True,270.28799999999995,4,3,COC(=O)CCc1nccc2c3ccc(OC)cc3n(C)c12,[0. 0. 0. ... 0. 0. 0.],mono +CC#CC1=CC=C(C2=CC=C(C(=O)O)S2)S1,True,248.32799999999997,2,2,COC(=O)c1ccc(-c2ccc(C3C=CC(C)CC3C)s2)s1,[0. 0. 0. ... 0. 0. 0.],mono +O=C1NCCC2=C1NC1=CC=CC=C21,True,186.21400000000003,0,2,Cn1c2c(c3ccccc31)CCNC2=O,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)C=CC1=CC=C(O)C2=C1C=C(C1=CC=C(O)C(O)=C1)O2,True,312.27700000000004,2,3,COC(=O)C=Cc1ccc(O)c2oc(-c3ccc(O)c(O)c3)cc12,[0. 0. 0. ... 0. 0. 0.],mono +CC#CC1=CC=C(C2=CC=C(C(=O)CCl)S2)S1,True,280.80099999999993,3,2,CCC(=O)c1ccc(-c2ccc(C3CC(C)C=CC3C)s2)s1,[0. 0. 0. ... 0. 0. 0.],mono +CC#CC1=CC=C(C2=CC=C(Cl)S2)S1,True,238.76399999999998,1,2,Cc1ccc(-c2ccc(C3CC(C)C=CC3C)s2)s1,[0. 0. 0. ... 0. 0. 0.],mono +O=C(O)[C@H]1O[C@@H](OC2=CC(O)=C3C(=O)C=C(C4=CC=CC=C4)OC3=C2)[C@H](O)[C@@H](O)[C@H]1O,True,430.36500000000007,3,3,COC(=O)[C@H]1O[C@@H](Oc2cc(O)c3c(=O)cc(-c4ccccc4)oc3c2)[C@H](O)[C@@H](O)[C@H]1O,[0. 0. 0. ... 0. 0. 0.],mono +N[C@@H](CC(=O)C1=CC=CC2=C1N=C1C(=CC(=O)C3=NC(C(=O)O)=CC(O)=C13)O2)C(=O)O,True,423.33700000000016,3,2,COC(=O)c1cc(O)c2c3nc4c(C(=O)C[C@H](N)C(=O)OC)cccc4oc-3cc(=O)c2n1,[0. 1. 0. ... 0. 0. 0.],mono +COC1=CC=C(C2=CC(=O)C3=C(O)C=C(O)C=C3O2)C=C1,True,284.267,2,3,COc1ccc(-c2cc(=O)c3c(O)cc(O)cc3o2)cc1,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC2=CC=C3C=C(O)C=CC3=C2C(OC)=C1O,True,270.28399999999993,2,3,COc1cc2ccc3cc(O)ccc3c2c(OC)c1O,[0. 0. 0. ... 0. 0. 0.],di +O=C1C=C(C2=CC=CC=C2)OC2=CC(O)=CC(O)=C12,True,254.24099999999999,1,3,O=c1cc(-c2ccccc2)oc2cc(O)cc(O)c12,[0. 0. 0. ... 0. 0. 0.],di +CC1=CC(O)=CC(O)=C1,True,124.13899999999997,0,1,Cc1cc(O)cc(O)c1,[0. 0. 0. ... 0. 0. 0.],di +O=C1C=CC2=CC(O)=C(O)C=C2O1,True,178.143,0,2,O=c1ccc2cc(O)c(O)cc2o1,[0. 0. 0. ... 0. 0. 0.],di +O=C1OC2=CC(O)=CC=C2C2=C1C1=CC=C(O)C=C1O2,True,268.224,0,4,O=c1oc2cc(O)ccc2c2oc3cc(O)ccc3c12,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=C(C2=COC3=CC(O)=CC(O)=C3C2=O)C=C1,True,284.267,2,3,COc1ccc(-c2coc3cc(O)cc(O)c3c2=O)cc1,[0. 0. 0. ... 0. 0. 0.],di +NC1=NC=NC2=C1N=CN2,True,135.13,0,2,Nc1ncnc2[nH]cnc12,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC(O)=C2C(=O)OC3=CC(O)=CC(C)=C3C2=C1,True,272.256,1,3,COc1cc(O)c2c(=O)oc3cc(O)cc(C)c3c2c1,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC(O)=C2C(O)=C3C(=O)C=C(C)OC3=CC2=C1,True,272.256,1,3,COc1cc(O)c2c(O)c3c(=O)cc(C)oc3cc2c1,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC(O)=CC2=CC=C3C=C(O)C=CC3=C12,True,240.25799999999998,1,3,COc1cc(O)cc2ccc3cc(O)ccc3c12,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC=CC(O)=C1,True,110.11199999999998,0,1,Oc1cccc(O)c1,[0. 0. 0. ... 0. 0. 0.],di +O=C1C2=CC=C(O)C=C2OC2=CC=CC(O)=C12,True,228.20299999999997,0,3,O=c1c2ccc(O)cc2oc2cccc(O)c12,[0. 0. 0. ... 0. 0. 0.],di +COC1=C(O)C=CC2=C1OC1=CC=CC(O)=C1C2=O,True,258.229,1,3,COc1c(O)ccc2c(=O)c3c(O)cccc3oc12,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC2=C(C=C1O)OC1=CC=CC(O)=C1C2=O,True,258.229,1,3,COc1cc2c(=O)c3c(O)cccc3oc2cc1O,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=C(O)C2=C1OC1=CC=C(O)C=C1C2=O,True,258.229,1,3,COc1ccc(O)c2c(=O)c3cc(O)ccc3oc12,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC=CC=C1O,True,110.11199999999998,0,1,Oc1ccccc1O,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC(O)=CC(C=CC2=CC=CC=C2)=C1,True,212.248,2,2,Oc1cc(O)cc(C=Cc2ccccc2)c1,[0. 0. 0. ... 0. 0. 0.],di +O=C1OC2=CC(O)=CC=C2C2=CC=C(O)C=C12,True,228.20299999999997,0,3,O=c1oc2cc(O)ccc2c2ccc(O)cc12,[0. 0. 0. ... 0. 0. 0.],di +COC1=C(O)C=C(O)C2=C1OC(C1=CC=CC=C1)=CC2=O,True,284.267,2,3,COc1c(O)cc(O)c2c(=O)cc(-c3ccccc3)oc12,[0. 0. 0. ... 0. 0. 0.],di +COC1=C(O)C=C2OC(C3=CC=CC=C3)=CC(=O)C2=C1O,True,284.267,2,3,COc1c(O)cc2oc(-c3ccccc3)cc(=O)c2c1O,[0. 0. 0. ... 0. 0. 0.],di +O=C1N=CNC2=C1NC=N2,True,136.11399999999998,0,2,O=c1nc[nH]c2nc[nH]c12,[0. 0. 0. ... 0. 0. 0.],di +CC1=CC2=NC3=C(N=C2C=C1C)C(=O)NC(=O)N3,True,242.23800000000003,0,3,Cc1cc2nc3[nH]c(=O)[nH]c(=O)c3nc2cc1C,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC2=C(C(O)=C1O)C(=O)C=C(C1=CC=CC=C1)O2,True,284.267,2,3,COc1cc2oc(-c3ccccc3)cc(=O)c2c(O)c1O,[0. 0. 0. ... 0. 0. 0.],di +O=C1C(C2=CC=C3OCOC3=C2)=COC2=CC(O)=CC(O)=C12,True,298.25,1,3,O=c1c(-c2ccc3c(c2)OCO3)coc2cc(O)cc(O)c12,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=C2OC3=CC(O)=CC(O)=C3C(=O)C2=C1,True,258.229,1,3,COc1ccc2oc3cc(O)cc(O)c3c(=O)c2c1,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC(O)=C2C(=O)C(O)=C(C3=CC=CC=C3)OC2=C1,True,284.26700000000005,2,3,COc1cc(O)c2c(=O)c(O)c(-c3ccccc3)oc2c1,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC(O)=C2C(=O)C3=C(C)C=C(O)C=C3OC2=C1,True,272.256,1,3,COc1cc(O)c2c(=O)c3c(C)cc(O)cc3oc2c1,[0. 0. 0. ... 0. 0. 0.],di +NC1=CC(Cl)=CC(Cl)=N1,True,163.007,0,1,Cc1cc(C)nc(N)c1,[0. 0. 0. ... 0. 0. 0.],di +O=C1C=CC2=C(O)C=C(O)C=C2O1,True,178.143,0,2,O=c1ccc2c(O)cc(O)cc2o1,[0. 0. 0. ... 0. 0. 0.],di +CC1=C(O)C=C(CCC2=CC=CC=C2)C=C1O,True,228.29100000000003,3,2,Cc1c(O)cc(CCc2ccccc2)cc1O,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=C2C=CC(=O)OC2=C1OC1=C(O)C=C(O)C2=C1OC(=O)C=C2,True,368.2970000000001,3,4,COc1ccc2ccc(=O)oc2c1Oc1c(O)cc(O)c2ccc(=O)oc12,[0. 0. 0. ... 0. 0. 0.],di +O=C1C=C(C2=CC=CC=C2)C2=CC(O)=C(O)C=C2O1,True,254.24099999999999,1,3,O=c1cc(-c2ccccc2)c2cc(O)c(O)cc2o1,[0. 0. 0. ... 0. 0. 0.],di +O=C1C=CC2=CC=C(O)C(O)=C2O1,True,178.14299999999997,0,2,O=c1ccc2ccc(O)c(O)c2o1,[0. 0. 0. ... 0. 0. 0.],di +ClC1=CC=C(Cl)C=C1,True,147.00400000000002,0,1,Cc1ccc(C)cc1,[0. 0. 0. ... 0. 0. 0.],di +ClC1=CC=CC=C1Cl,True,147.00399999999996,0,1,Cc1ccccc1C,[0. 0. 0. ... 0. 0. 0.],di +CC1=CNN=C1,True,82.106,0,1,Cc1cn[nH]c1,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC(O)=CC(C2=CC3=CC4=C(C=C3O2)OC=C4)=C1,True,266.25199999999995,1,4,Oc1cc(O)cc(-c2cc3cc4ccoc4cc3o2)c1,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC(O)=C2C(=O)C3=CC=CC(O)=C3OC2=C1,True,258.229,1,3,COc1cc(O)c2c(=O)c3cccc(O)c3oc2c1,[0. 0. 0. ... 0. 0. 0.],di +C=CC1=CC(O)=CC2=CC=C3C(C)=C(O)C=CC3=C12,True,250.29699999999997,1,3,C=Cc1cc(O)cc2ccc3c(C)c(O)ccc3c12,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC(O)=C2C=CC=CC2=C1,True,160.17199999999997,0,2,Oc1cc(O)c2ccccc2c1,[0. 0. 0. ... 0. 0. 0.],di +O=C1C2=CC(O)=CC=C2OC2=CC3=C(OCO3)C(O)=C12,True,272.21199999999993,0,3,O=c1c2cc(O)ccc2oc2cc3c(c(O)c12)OCO3,[0. 0. 0. ... 0. 0. 0.],di +COC1=C(O)C=C2C=CC3=CC(O)=CC=C3C2=C1OC,True,270.284,2,3,COc1c(O)cc2ccc3cc(O)ccc3c2c1OC,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC(O)=CC2=CC=C3C=CC=C(O)C3=C12,True,240.25799999999998,1,3,COc1cc(O)cc2ccc3cccc(O)c3c12,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=C(C2=COC3=CC(O)=C(O)C=C3C2=O)C=C1,True,284.267,2,3,COc1ccc(-c2coc3cc(O)c(O)cc3c2=O)cc1,[0. 0. 0. ... 0. 0. 0.],di +O=C1C=C(C2=CC=CC=C2)OC2=C(O)C=CC(O)=C12,True,254.24099999999999,1,3,O=c1cc(-c2ccccc2)oc2c(O)ccc(O)c12,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=C(OC2=CC(=O)C3=C(O)C=C(O)C=C3O2)C=C1,True,300.266,3,3,COc1ccc(Oc2cc(=O)c3c(O)cc(O)cc3o2)cc1,[0. 0. 0. ... 0. 0. 0.],di +O=C1C2=CC=CC=C2C(=O)C2=C1C=C(O)C(CC1=CC=CC=C1)=C2O,True,330.339,2,3,O=C1c2ccccc2C(=O)c2c1cc(O)c(Cc1ccccc1)c2O,[0. 1. 0. ... 0. 0. 0.],di +OC1=CC=C(O)C=C1,True,110.11199999999998,0,1,Oc1ccc(O)cc1,[0. 0. 0. ... 0. 0. 0.],di +CC1=CC=C(O)C(O)=C1,True,124.13899999999998,0,1,Cc1ccc(O)c(O)c1,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=CC=C1C1=COC2=CC(O)=C(C)C(O)=C2C1=O,True,298.29400000000004,2,3,COc1ccccc1-c1coc2cc(O)c(C)c(O)c2c1=O,[0. 0. 0. ... 0. 0. 0.],di +C=CC1=C(C)C(O)=CC2=CC=C3C(C)=C(O)C=CC3=C12,True,264.324,1,3,C=Cc1c(C)c(O)cc2ccc3c(C)c(O)ccc3c12,[0. 0. 0. ... 0. 0. 0.],di +O=C1C=C(CCC2=CC=CC=C2)OC2=C(O)C=CC(O)=C12,True,282.295,3,3,O=c1cc(CCc2ccccc2)oc2c(O)ccc(O)c12,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC=C(O)C2=CC=CC=C12,True,160.172,0,2,Oc1ccc(O)c2ccccc12,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC(O)=CC(CCC2=CC=CC=C2)=C1,True,214.264,3,2,Oc1cc(O)cc(CCc2ccccc2)c1,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC(O)=C2OC3=CC=CC=C3CCC2=C1,True,228.24699999999996,0,2,Oc1cc(O)c2c(c1)CCc1ccccc1O2,[0. 0. 0. ... 0. 0. 0.],di +CC1=CC(O)=CC=C1O,True,124.13899999999997,0,1,Cc1cc(O)ccc1O,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC(O)=C2C(=C1)C=CC1=CC(O)=CC=C12,True,240.258,1,3,COc1cc(O)c2c(ccc3cc(O)ccc32)c1,[0. 0. 0. ... 0. 0. 0.],di +COC1=C(O)C=C2OC3=CC=CC=C3C(=O)C2=C1O,True,258.229,1,3,COc1c(O)cc2oc3ccccc3c(=O)c2c1O,[0. 0. 0. ... 0. 0. 0.],di +CC1=CC=CC(O)=C1O,True,124.13899999999998,0,1,Cc1cccc(O)c1O,[0. 0. 0. ... 0. 0. 0.],di +O=C1OC2=CC(O)=C(O)C=C2C2=C1C1=CC=CC=C1O2,True,268.224,0,4,O=c1oc2cc(O)c(O)cc2c2oc3ccccc3c12,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC2=CC=C3C=C(O)C=CC3=C2C(O)=C1OC,True,270.284,2,3,COc1cc2ccc3cc(O)ccc3c2c(O)c1OC,[0. 0. 0. ... 0. 0. 0.],di +NC1=NC=NC2=C1NC=N2,True,135.12999999999997,0,2,Nc1ncnc2nc[nH]c12,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=C(C2=COC3=C(O)C(O)=CC=C3C2=O)C=C1,True,284.267,2,3,COc1ccc(-c2coc3c(O)c(O)ccc3c2=O)cc1,[0. 0. 0. ... 0. 0. 0.],di +CC1=CC(=O)OC2=CC(O)=C(O)C=C12,True,192.17,0,2,Cc1cc(=O)oc2cc(O)c(O)cc12,[0. 0. 0. ... 0. 0. 0.],di +CC1=C(O)C=C(C=CC2=CC=CC=C2)C=C1O,True,226.275,2,2,Cc1c(O)cc(C=Cc2ccccc2)cc1O,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC=C2C=CC=CC2=C1O,True,160.17199999999997,0,2,Oc1ccc2ccccc2c1O,[0. 0. 0. ... 0. 0. 0.],di +OC1=CN=C2C(O)=CC=CC2=C1,True,161.15999999999997,0,2,Oc1cnc2c(O)cccc2c1,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC2=C(C=C1O)NC=C2,True,149.149,0,2,Cn1ccc2cc(O)c(O)cc21,[0. 0. 0. ... 0. 0. 0.],di +O=C1C(C2=CC=CC=C2)=COC2=CC(O)=CC(O)=C12,True,254.24099999999996,1,3,O=c1c(-c2ccccc2)coc2cc(O)cc(O)c12,[0. 0. 0. ... 0. 0. 0.],di +O=C1C=C(C2=CC=CC=C2)C2=CC=C(O)C(O)=C2O1,True,254.24099999999999,1,3,O=c1cc(-c2ccccc2)c2ccc(O)c(O)c2o1,[0. 0. 0. ... 0. 0. 0.],di +CC1=C(O)C=CC2=C1C=CC1=CC(O)=CC(C=O)=C12,True,252.26899999999998,1,3,Cc1c(O)ccc2c1ccc1cc(O)cc(C=O)c12,[0. 0. 0. ... 0. 0. 0.],di +COC1=C(O)C=CC2=C1OC(=O)C1=C2OC2=CC(O)=CC=C21,True,298.24999999999994,1,4,COc1c(O)ccc2c1oc(=O)c1c3ccc(O)cc3oc21,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=C2C(=C1)OCC1=C2OC2=CC(O)=C(O)C=C12,True,284.267,1,3,COc1ccc2c(c1)OCc1c-2oc2cc(O)c(O)cc12,[0. 0. 0. ... 0. 0. 0.],di +CC1=C(O)C=C(C2=CC3=CC=CC=C3O2)C=C1O,True,240.258,1,3,Cc1c(O)cc(-c2cc3ccccc3o2)cc1O,[0. 0. 0. ... 0. 0. 0.],di +CC1=C(O)C(C=O)=C2OC(C3=CC=CC=C3)=CC(=O)C2=C1O,True,296.278,2,3,Cc1c(O)c(C=O)c2oc(-c3ccccc3)cc(=O)c2c1O,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=C(O)C2=C1OC1=CC(C)=CC(O)=C1C2=O,True,272.25600000000003,1,3,COc1ccc(O)c2c(=O)c3c(O)cc(C)cc3oc12,[0. 0. 0. ... 0. 0. 0.],di +COC1=C(C)C(O)=C2C(=O)C=C(C3=CC=CC=C3)OC2=C1O,True,298.29400000000004,2,3,COc1c(C)c(O)c2c(=O)cc(-c3ccccc3)oc2c1O,[0. 0. 0. ... 0. 0. 0.],di +COC1=C(O)C=C2C(=O)OC3=CC(O)=CC4=CC=C1C2=C43,True,282.25100000000003,1,4,COc1c(O)cc2c(=O)oc3cc(O)cc4ccc1c2c43,[0. 0. 0. ... 0. 0. 0.],di +CC1=CC2=C(C(=O)O1)C1=CC(O)=C(O)C=C1C(=O)O2,True,260.201,0,3,Cc1cc2oc(=O)c3cc(O)c(O)cc3c2c(=O)o1,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC=C2C(O)=CC=CC2=C1,True,160.172,0,2,Oc1ccc2c(O)cccc2c1,[0. 0. 0. ... 0. 0. 0.],di +O=C1C=C(CCC2=CC=CC=C2)OC2=C(O)C=C(O)C=C12,True,282.295,3,3,O=c1cc(CCc2ccccc2)oc2c(O)cc(O)cc12,[0. 0. 0. ... 0. 0. 0.],di +COC1=CC=C2C(=C1)OC1=C2C(=O)OC2=CC(O)=C(O)C=C21,True,298.25,1,4,COc1ccc2c(c1)oc1c3cc(O)c(O)cc3oc(=O)c21,[0. 0. 0. ... 0. 0. 0.],di +NC1=NC=C2N=C(O)C=NC2=N1,True,163.14,0,2,Nc1ncc2nc(O)cnc2n1,[0. 0. 0. ... 0. 0. 0.],di +COC1=C(O)C=CC2=C1C(=O)C1=CC(O)=CC=C1O2,True,258.229,1,3,COc1c(O)ccc2oc3ccc(O)cc3c(=O)c12,[0. 0. 0. ... 0. 0. 0.],di +OC1=CC=C2C=C(O)C=CC2=C1,True,160.172,0,2,Oc1ccc2cc(O)ccc2c1,[0. 0. 0. ... 0. 0. 0.],di diff --git a/Examples/example_calculation_submission.py b/Examples/example_g16_calculation_submission.py similarity index 100% rename from Examples/example_calculation_submission.py rename to Examples/example_g16_calculation_submission.py diff --git a/Examples/example_psi4_opt_submission.py b/Examples/example_psi4_opt_submission.py new file mode 100644 index 0000000..0c0d5bf --- /dev/null +++ b/Examples/example_psi4_opt_submission.py @@ -0,0 +1,21 @@ +import rdkit +from rdkit import Chem +from rdkit.Chem import Draw +from rdkit.Chem import AllChem +import numpy as np +import pandas as pd +import QCflow +from QCflow.load_gaussian import * +from QCflow.torsion_parser import * +from QCflow.find_torsion import * +from QCflow.write_psi4 import * +from QCflow.run_psi4 import * +from QCflow.energy_calculations import * +import os + +df = pd.read_csv('FilteredSmilesFragments_props.csv') + +smiles_dict = {str(i): row['smiles'] for i, row in df.iterrows()} + +for k, v in smiles_dict.items(): + run_psi4('opt', k, v, 12, 10, 'b3lyp', '3-21G') diff --git a/Examples/example_torsional_scan_script.py b/Examples/example_torsional_scan_script.py new file mode 100644 index 0000000..7872cc4 --- /dev/null +++ b/Examples/example_torsional_scan_script.py @@ -0,0 +1,43 @@ +import rdkit +from rdkit import Chem +from rdkit.Chem import Draw +from rdkit.Chem import AllChem +import numpy as np +import pandas as pd +from QCflow.load_gaussian import * +from QCflow.torsion_parser import * +from QCflow.find_torsion import * +from QCflow.write_psi4 import * +from QCflow.run_psi4 import * +from QCflow.energy_calculations import * +import os + +# This is an example for one molecule +# This can easily be converted into a loop for multiple molecules +# The getTorsion assumes single bonded but can be altered if required +# If you have more than 1 torsion or tiple, double etc bond link +# Refer to this --> https://github.com/matta-research-group/QCflow/blob/qcflow-0.5/QCflow/torsion_parser.py +# torsion_parser.py also has the scripts to extract the data you may want +# The docs --> https://matta-research-group.github.io/QCflow/ + +smiles_example = 'CC(=O)Oc1c(OC(C)=O)c(-c2c(OC(C)=O)c(OC(C)=O)cc3cc[nH]c23)c2cc[nH]c2c1' + +mol = Chem.MolFromSmiles(smiles_example) + +#adds H's and 3D coordinates +mol = embed_molecule(mol) + +#finds the rotatable bond in the molecule +bond = getBond(mol) + +#Gets the whole torsion +mol_torsions = getTorsion(mol, bond[0]) + +#writes the Gaussian input file for the torsion scan +write_gaussian('tor', 'example_mol', smiles_example, functional='B3LYP', basis_set='6-31G*', mol=mol, + torsion=mol_torsions) +#writes the SLURM script to run the Gaussian job +write_slurm('tor', 'example_mol', cpus=10) + +#submits the SLURM job, with retries in case of failure +submit_slurm_job('tor', 'example_mol', max_retries=5, wait_seconds=30) \ No newline at end of file diff --git a/QCflow.yml b/QCflow.yml index e06a608..af1fa1e 100644 --- a/QCflow.yml +++ b/QCflow.yml @@ -1,9 +1,7 @@ name: QCflow channels: - conda-forge - - lich - - defaults - - anaconda + # - lich causing issues in instal dependencies: - brotli - brotli-bin @@ -103,7 +101,6 @@ dependencies: - seaborn-base - six - statsmodels - - syba - tk - tornado - typing-extensions diff --git a/QCflow/__init__.py b/QCflow/__init__.py index 72ea573..77dbda4 100755 --- a/QCflow/__init__.py +++ b/QCflow/__init__.py @@ -10,3 +10,5 @@ from . import write_psi4 from . import run_psi4 from . import sa_score +from . import write_xTB +from . import run_xTB diff --git a/QCflow/__pycache__/__init__.cpython-313.pyc b/QCflow/__pycache__/__init__.cpython-313.pyc index 6d0076c..c60f286 100644 Binary files a/QCflow/__pycache__/__init__.cpython-313.pyc and b/QCflow/__pycache__/__init__.cpython-313.pyc differ diff --git a/QCflow/__pycache__/energy_calculations.cpython-313.pyc b/QCflow/__pycache__/energy_calculations.cpython-313.pyc index 94d239d..fcaa04d 100644 Binary files a/QCflow/__pycache__/energy_calculations.cpython-313.pyc and b/QCflow/__pycache__/energy_calculations.cpython-313.pyc differ diff --git a/QCflow/__pycache__/find_torsion.cpython-313.pyc b/QCflow/__pycache__/find_torsion.cpython-313.pyc index b26e99f..b2393fe 100644 Binary files a/QCflow/__pycache__/find_torsion.cpython-313.pyc and b/QCflow/__pycache__/find_torsion.cpython-313.pyc differ diff --git a/QCflow/__pycache__/fragments.cpython-313.pyc b/QCflow/__pycache__/fragments.cpython-313.pyc index 6c8da25..ef10061 100644 Binary files a/QCflow/__pycache__/fragments.cpython-313.pyc and b/QCflow/__pycache__/fragments.cpython-313.pyc differ diff --git a/QCflow/__pycache__/load_gaussian.cpython-313.pyc b/QCflow/__pycache__/load_gaussian.cpython-313.pyc index 2e81a93..45bc0dc 100644 Binary files a/QCflow/__pycache__/load_gaussian.cpython-313.pyc and b/QCflow/__pycache__/load_gaussian.cpython-313.pyc differ diff --git a/QCflow/__pycache__/run_gaussian.cpython-313.pyc b/QCflow/__pycache__/run_gaussian.cpython-313.pyc index 25cd30d..80f345e 100644 Binary files a/QCflow/__pycache__/run_gaussian.cpython-313.pyc and b/QCflow/__pycache__/run_gaussian.cpython-313.pyc differ diff --git a/QCflow/__pycache__/run_psi4.cpython-313.pyc b/QCflow/__pycache__/run_psi4.cpython-313.pyc index e1a8511..8df41d4 100644 Binary files a/QCflow/__pycache__/run_psi4.cpython-313.pyc and b/QCflow/__pycache__/run_psi4.cpython-313.pyc differ diff --git a/QCflow/__pycache__/run_xTB.cpython-313.pyc b/QCflow/__pycache__/run_xTB.cpython-313.pyc new file mode 100644 index 0000000..5b83dd6 Binary files /dev/null and b/QCflow/__pycache__/run_xTB.cpython-313.pyc differ diff --git a/QCflow/__pycache__/sa_score.cpython-313.pyc b/QCflow/__pycache__/sa_score.cpython-313.pyc index 5c48f6e..83e41af 100644 Binary files a/QCflow/__pycache__/sa_score.cpython-313.pyc and b/QCflow/__pycache__/sa_score.cpython-313.pyc differ diff --git a/QCflow/__pycache__/slurm.cpython-313.pyc b/QCflow/__pycache__/slurm.cpython-313.pyc index 8f7abba..fe4506b 100644 Binary files a/QCflow/__pycache__/slurm.cpython-313.pyc and b/QCflow/__pycache__/slurm.cpython-313.pyc differ diff --git a/QCflow/__pycache__/testing_data.cpython-313.pyc b/QCflow/__pycache__/testing_data.cpython-313.pyc index f4e3baf..709669e 100644 Binary files a/QCflow/__pycache__/testing_data.cpython-313.pyc and b/QCflow/__pycache__/testing_data.cpython-313.pyc differ diff --git a/QCflow/__pycache__/torsion_parser.cpython-313.pyc b/QCflow/__pycache__/torsion_parser.cpython-313.pyc index 5ff8b48..f9fbf22 100644 Binary files a/QCflow/__pycache__/torsion_parser.cpython-313.pyc and b/QCflow/__pycache__/torsion_parser.cpython-313.pyc differ diff --git a/QCflow/__pycache__/write_gaussian.cpython-313.pyc b/QCflow/__pycache__/write_gaussian.cpython-313.pyc index eaf4f7a..ef76f1c 100644 Binary files a/QCflow/__pycache__/write_gaussian.cpython-313.pyc and b/QCflow/__pycache__/write_gaussian.cpython-313.pyc differ diff --git a/QCflow/__pycache__/write_psi4.cpython-313.pyc b/QCflow/__pycache__/write_psi4.cpython-313.pyc index c8a78d5..2762c1e 100644 Binary files a/QCflow/__pycache__/write_psi4.cpython-313.pyc and b/QCflow/__pycache__/write_psi4.cpython-313.pyc differ diff --git a/QCflow/__pycache__/write_xTB.cpython-313.pyc b/QCflow/__pycache__/write_xTB.cpython-313.pyc new file mode 100644 index 0000000..1f180e7 Binary files /dev/null and b/QCflow/__pycache__/write_xTB.cpython-313.pyc differ diff --git a/QCflow/run_psi4.py b/QCflow/run_psi4.py index 417138f..7009f82 100644 --- a/QCflow/run_psi4.py +++ b/QCflow/run_psi4.py @@ -17,6 +17,7 @@ def run_psi4(job_name, mol_name, mol_smile, time=24, cpus=10, functional='b3lyp' job_name (str): The type of job to run. Possible values: - 'sp': Single Point neutral - 'opt': Optimisation neutral + - 'opt_pre_geom': Optimisation neutral where a txt file called '{mol_name}_opt.xyz' is present in the molecule directory. This file should contain the geometry to be used for the optimisation in XYZ format. If this file is not present, the function will default to 'opt'. - 'cation': Geometry optimisation cation (opt_c) and single of neutral charge, cation geometry (n_c_geo) - 'anion': Geometry optimisation anion (opt_a) and single of neutral charge, anion geometry (n_a_geo) - 'sp_c': Single point calculation of neutral geometry at cation charge @@ -50,7 +51,14 @@ def run_psi4(job_name, mol_name, mol_smile, time=24, cpus=10, functional='b3lyp' os.mkdir(f'{mol_name}') #makes a directory for the molecule os.chdir(f'{mol_name}') #goes into directory - if (job_name=='opt') or (job_name=='sp'): + if (job_name=='opt') or (job_name=='sp') or (job_name=='opt_pre_geom'): + + if (job_name=='opt_pre_geom'): + if os.path.exists(f'{mol_name}_opt.xyz'): + print(f"Found {mol_name}_opt.xyz, using this geometry for optimisation.") + else: + print(f"{mol_name}_opt.xyz not found, defaulting to 'opt' job type.") + job_name = 'opt' #turns smiles string into rdkit object mol = Chem.MolFromSmiles(mol_smile) #gets rdkit estimated coordinates of dimer @@ -64,6 +72,9 @@ def run_psi4(job_name, mol_name, mol_smile, time=24, cpus=10, functional='b3lyp' if (job_name=='cation') or (job_name=='anion') or (job_name=='sp_c') or (job_name=='sp_a'): write_psi4_reorg(job_name, mol_name, functional, basis_set) + + else: + print(f"Invalid job_name: {job_name}, please use 'opt', 'sp', 'opt_pre_geom', 'cation', 'anion', 'sp_c' or 'sp_a'. See documentation for more details.") #writes the slurm file diff --git a/QCflow/run_xTB.py b/QCflow/run_xTB.py new file mode 100644 index 0000000..0a88d28 --- /dev/null +++ b/QCflow/run_xTB.py @@ -0,0 +1,76 @@ +from QCflow.fragments import * +from QCflow.find_torsion import * +from QCflow.write_gaussian import * +from QCflow.slurm import * +from QCflow.torsion_parser import * +from QCflow.run_gaussian import * +from QCflow.write_psi4 import * +from QCflow.write_xTB import * +import os +import json +import subprocess +import tempfile + +def run_xTB(job_name, mol_name, mol_smile, time=2, cpus=10, functional='b3lyp', basis_set='6-31g*'): + """ + Submits a xTB calculation to CREATE HPC. + + Parameters + ---------- + job_name (str): The type of job to run. Possible values: + - 'opt': Optimisation neutral + - 'cation': Geometry optimisation cation (opt_c) and single of neutral charge, cation geometry (n_c_geo) + - 'anion': Geometry optimisation anion (opt_a) and single of neutral charge, anion geometry (n_a_geo) + - 'sp_c': Single point calculation of neutral geometry at cation charge + - 'sp_a': Single point calculation of neutral geometry at anion charge + mol_name : str + Name of the oligomer. + mol_smile : str + SMILES string of the oligomer. + time (int, optional): The time limit for the job in hours. Default is 24. (Max is 48) + cpus (int, optional): The number of CPUs to allocate for the job. Default is 10. + functional : str, optional + Quantum chemistry functional to be used (default is 'B3LYP'). + basis_set : str, optional + Basis set to be used (default is '6-31G*'). + + + Returns + ------- + - For 'opt': + - Converts the SMILES string to an RDKit molecule object. + - Generates 3D coordinates for the molecule. + - Predicts the conformer geometry. + - Runs xTB optimisation using the generated coordinates. + - Writes a psi4 input file with the optimised geometry. + - Runs a single point calculation. + Finally, the function writes a SLURM script, submits the job, and returns to the previous directory. + """ + if os.path.exists(f'{mol_name}'): + os.chdir(f'{mol_name}') #goes into directory + else: + os.mkdir(f'{mol_name}') #makes a directory for the molecule + os.chdir(f'{mol_name}') #goes into directory + + if (job_name=='opt'): + #turns smiles string into rdkit object + mol = Chem.MolFromSmiles(mol_smile) + #gets rdkit estimated coordinates of dimer + mol3d = embed_molecule(mol) + + conf_geo = rdkit_predict_conf(mol_smile) + #writes a guassian input file + write_xTB_psi4(job_name, mol_name, mol_smile, functional, basis_set, mol=mol3d, conformer=conf_geo) + + #for reorganisation calcultions + if (job_name=='cation') or (job_name=='anion') or (job_name=='sp_c') or (job_name=='sp_a'): + + write_xTB_psi4_reorg(job_name, mol_name, functional, basis_set) + + + #writes the slurm file + write_slurm_psi4(job_name, mol_name, time, cpus) + #submits the slurm jon + submit_slurm_job(job_name, mol_name) + #goes back to previous dirctory + os.chdir(os.path.dirname(os.getcwd())) diff --git a/QCflow/write_psi4.py b/QCflow/write_psi4.py index 322f9b0..92c8746 100755 --- a/QCflow/write_psi4.py +++ b/QCflow/write_psi4.py @@ -20,7 +20,8 @@ def write_psi4(job_name, mol_name, smile, functional='b3lyp', basis_set='6-31g*' ---------- job_name (str): The type of job to run. Possible values: - 'sp': Single Point neutral - - 'opt': Optimisation neutral + - 'opt': Optimisation neutral + - 'opt_pre_geom' : Optimisation neutral, from stated geometry, where a txt file called '{mol_name}_opt.xyz' is present and in psi4 format mol_name (str): The name of the molecule. smile (str): The SMILE string of the molecule. functional (str, optional): The functional to use. Default is 'b3lyp'. @@ -66,7 +67,29 @@ def write_psi4(job_name, mol_name, smile, functional='b3lyp', basis_set='6-31g*' if Descriptors.NumRadicalElectrons(mol) == 1: mult_chg = '0 2' if Descriptors.NumRadicalElectrons(mol) == 2: - mult_chg = '0 3' + mult_chg = '0 3' + + if (job_name=='opt_pre_geom'): + + old_chk = ' ' + + geometry_file = f'{mol_name}_opt.xyz' + calculation = 'opt' + torsion_data = f' \n' + + if Descriptors.NumRadicalElectrons(mol) == 0: + charge = '0' + mult = '1' + mult_chg = f'{charge} {mult}' + if Descriptors.NumRadicalElectrons(mol) == 1: + charge = '0' + mult = '2' + mult_chg = f'{charge} {mult}' + if Descriptors.NumRadicalElectrons(mol) == 2: + charge = '0' + mult = '3' + mult_chg = f'{charge} {mult}' + file_name = f'{mol_name}_{job_name}.py' #Includes information about basis set to allow for ramping @@ -74,56 +97,98 @@ def write_psi4(job_name, mol_name, smile, functional='b3lyp', basis_set='6-31g*' number_proc = 'set_num_threads(4)' memory_num = '4000 mb' - - with open(file_name, 'w') as file: - file.write('import psi4\n') - file.write(' \n')# - file.write(f"psi4.set_memory('{memory_num}')\n") - file.write('psi4.set_num_threads(4)\n') - file.write(' \n')# - file.write('mol_name = psi4.geometry("""\n')# - #file.write(' symmetry c1 \n') - file.write(f' {mult_chg} \n')# - # if no previous checkpoint file, get geometry - if old_chk==' ': - for atom,symbol in enumerate(symbols): - p = geometry.GetAtomPosition(atom) - # atom x y z - line = f' {symbol} {p.x:.5f} {p.y:.5f} {p.z:.5f} \n' - file.write(line) - file.write('""") \n') - file.write(' \n') - file.write("psi4.set_module_options('scf', {'maxiter': 75}) \n") - file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df'}})\n") - file.write(' \n') - if (job_name=='sp'): - file.write(f"energy, wfn = psi4.energy('{functional}', return_wfn=True) \n") - file.write(' \n') - file.write("sp_geometry_xyz = wfn.molecule().to_string('xyz') \n") - file.write(f"with open('{mol_name}_{job_name}.xyz', 'w') as xyz_file:\n") - file.write(" xyz_file.write(sp_geometry_xyz) \n") - file.write(' \n') - file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") - file.write(' \n') - file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') - file.write(' \n') - file.write('n_occ_alpha = wfn.nalpha() \n') - file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') - file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') - file.write('energy_gap = lumo_energy - homo_energy \n') + if (job_name== 'sp') or (job_name=='opt'): + with open(file_name, 'w') as file: + file.write('import psi4\n') + file.write(' \n')# + file.write(f"psi4.set_memory('{memory_num}')\n") + file.write('psi4.set_num_threads(4)\n') + file.write(' \n')# + file.write('mol_name = psi4.geometry("""\n')# + #file.write(' symmetry c1 \n') + file.write(f' {mult_chg} \n')# + # if no previous checkpoint file, get geometry + if old_chk==' ': + for atom,symbol in enumerate(symbols): + p = geometry.GetAtomPosition(atom) + # atom x y z + line = f' {symbol} {p.x:.5f} {p.y:.5f} {p.z:.5f} \n' + file.write(line) + file.write('""") \n') + file.write(' \n') + file.write("psi4.set_options({'SCF__MAXITER': 100, 'GEOM_MAXITER': 200}) \n") + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df'}})\n") file.write(' \n') - file.write('hartree_to_ev = 27.2114079527 \n') - file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') - file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') - file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + if (job_name=='sp'): + file.write(f"energy, wfn = psi4.energy('{functional}', return_wfn=True) \n") + file.write(' \n') + file.write("sp_geometry_xyz = wfn.molecule().to_string('xyz') \n") + file.write(f"with open('{mol_name}_{job_name}.xyz', 'w') as xyz_file:\n") + file.write(" xyz_file.write(sp_geometry_xyz) \n") + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('hartree_to_ev = 27.2114079527 \n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_{job_name}_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Optimized energy: {energy:.6f} Hatree\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + if (job_name=='opt'): + file.write(f"energy, wfn = psi4.optimize('{functional}', engine='geometric', return_wfn=True) \n") + file.write(' \n') + file.write("optimized_geometry_xyz = wfn.molecule().to_string('xyz') \n") + file.write(f"with open('{mol_name}_{job_name}.xyz', 'w') as xyz_file:\n") + file.write(" xyz_file.write(optimized_geometry_xyz) \n") + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('hartree_to_ev = 27.2114079527 \n') + file.write('energy_ev = energy * hartree_to_ev \n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_{job_name}_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Optimized energy: {energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + file.write(' \n') + + if (job_name == 'opt_pre_geom'): + with open(file_name, 'w') as file: + file.write('import psi4\n') + file.write(' \n')# + file.write(f"psi4.set_memory('{memory_num}')\n") + file.write('psi4.set_num_threads(4)\n') + file.write(' \n')# + file.write(f"mol = psi4.geometry(open('{geometry_file}').read())\n")#use geometry from opt + file.write('psi4.core.set_active_molecule(mol)\n')#set as active molecule + file.write(f"mol.set_molecular_charge({charge}) \n")# set charge + file.write(f"mol.set_multiplicity({mult}) \n")# set multiplicity + file.write("psi4.set_options({'SCF__MAXITER': 100, 'GEOM_MAXITER': 200}) \n") + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df'}})\n") file.write(' \n') - file.write(f"with open('{mol_name}_{job_name}_energy_and_gap.txt', 'w') as file:\n") - file.write(' file.write(f"Optimized energy: {energy:.6f} Hatree\\n") \n') - file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') - file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') - file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') - if (job_name=='opt'): - file.write(f"energy, wfn = psi4.optimize('{functional}', return_wfn=True) \n") + file.write(f"energy, wfn = psi4.optimize('{functional}', engine='geometric', return_wfn=True) \n") file.write(' \n') file.write("optimized_geometry_xyz = wfn.molecule().to_string('xyz') \n") file.write(f"with open('{mol_name}_{job_name}.xyz', 'w') as xyz_file:\n") @@ -149,7 +214,10 @@ def write_psi4(job_name, mol_name, smile, functional='b3lyp', basis_set='6-31g*' file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') - file.write(' \n') + file.write(' \n') + + else: + print(f"Invalid job_name: {job_name}, please use 'opt', 'sp' or 'opt_pre_geom'. See documentation for more details.") def write_psi4_reorg(job_name, mol_name, functional='b3lyp', basis_set='6-31g*'): @@ -227,10 +295,10 @@ def write_psi4_reorg(job_name, mol_name, functional='b3lyp', basis_set='6-31g*') file.write(f"mol.set_multiplicity({mult}) \n")# set multiplicity file.write(' \n') if (job_name=='cation') or (job_name=='anion'): - file.write("psi4.set_module_options('scf', {'maxiter': 100}) \n") + file.write("psi4.set_options({'SCF__MAXITER': 100, 'GEOM_MAXITER': 200}) \n") file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df', 'reference': 'uhf'}})\n") file.write(' \n') - file.write(f"energy, wfn = psi4.optimize('{functional}', return_wfn=True) \n") + file.write(f"energy, wfn = psi4.optimize('{functional}', engine='geometric', return_wfn=True) \n") file.write(' \n') file.write("optimized_geometry_xyz = wfn.molecule().to_string('xyz') \n") file.write(f"with open('{mol_name}_{job_type}.xyz', 'w') as xyz_file:\n") @@ -316,5 +384,7 @@ def write_psi4_reorg(job_name, mol_name, functional='b3lyp', basis_set='6-31g*') file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') file.write(' \n') + else: + print(f"Invalid job_name: {job_name}, please use 'cation', 'anion', 'sp_a' or 'sp_c'. See documentation for more details.") \ No newline at end of file diff --git a/QCflow/write_xTB.py b/QCflow/write_xTB.py new file mode 100644 index 0000000..e921a22 --- /dev/null +++ b/QCflow/write_xTB.py @@ -0,0 +1,376 @@ +from QCflow.fragments import * +from QCflow.find_torsion import * +from QCflow.write_gaussian import * +from QCflow.slurm import * +from QCflow.torsion_parser import * +from QCflow.run_gaussian import * +from QCflow.write_psi4 import * +import os +import json +import subprocess +import tempfile + +def run_xtb_optimisation(xyz: str, output_file: str, charge=0, uhf=0, gfn=2): + """ + Runs a xTB geometry optimisation using GFN2-xTB and saves the result to output_file. + + Parameters + ---------- + xyz (str): The input geometry in XYZ format as a string. + output_file (str): The path to the file where the optimised geometry will be saved + charge (int, optional): The total charge of the system. Default is 0. + uhf (int, optional): The number of unpaired electrons (spin multiplicity - 1). Default is 0 (singlet). + gfn (int, optional): The GFN-xTB method to use (1, 2, or 3). Default is 2 (GFN2-xTB). + + Returns + ------- + None: Writes the optimised geometry to output_file. + + """ + with tempfile.TemporaryDirectory() as tmpdir: + xyz_path = os.path.join(tmpdir, "input.xyz") + with open(xyz_path, 'w') as f: + f.write(xyz) + + cmd = [ + "xtb", xyz_path, + "--opt", + "--chrg", str(charge), + "--uhf", str(uhf), + "--gfn", str(gfn) + ] + + try: + subprocess.run(cmd, cwd=tmpdir, check=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) + except subprocess.CalledProcessError as e: + raise RuntimeError(f"xTB failed: {e.stderr.decode()}") + + opt_path = os.path.join(tmpdir, "xtbopt.xyz") + if not os.path.exists(opt_path): + raise FileNotFoundError("Optimised structure (xtbopt.xyz) not found.") + + with open(opt_path, 'r') as f: + optimised_xyz = f.read() + + with open(output_file, 'w') as f: + f.write(optimised_xyz) + + print(f"Optimised structure saved to: {output_file}") + +def load_xyz_from_file(filepath): + """ + Loads an XYZ file and returns the content as a string. + + Parameters + ---------- + filepath (str): The path to the XYZ file to be loaded. + + Returns + ------- + str: The content of the XYZ file as a string. + """ + with open(filepath, 'r') as f: + lines = f.readlines() + + # Remove the first two lines (atom count and comment/energy) + coordinates = lines[2:] + + # Recalculate number of atoms just in case + num_atoms = len(coordinates) + + # Build the xyz string + xyz_string = f"{num_atoms}\n\n" + "".join(coordinates) + return xyz_string + + +def write_xTB_psi4(job_name, mol_name, smile, functional='b3lyp', basis_set='6-31g*', mol=None, conformer=None): + """ + Runs a xTB geometry optimisation and caulcates electronic proeprities using Psi4, based on the provided parameters. + + Parameters + ---------- + job_name (str): The type of job to run. Possible values: + - 'opt': Optimisation neutral + mol_name (str): The name of the molecule. + smile (str): The SMILE string of the molecule. + functional (str, optional): The functional to use. Default is 'b3lyp'. + basis_set (str, optional): The basis set to use. Default is '6-31g'. + mol (rdkit.Chem.rdchem.Mol, optional): The RDKit embedded molecule object. + conformer (rdkit.Chem.rdchem.Conformer, optional): The RDKit conformer of the molecule. + + Returns + ------- + None: Writes the psi4 input file to disk. + + """ + if (job_name=='opt'): + + old_chk = ' ' + symbols = [a.GetSymbol() for a in mol.GetAtoms()] + # get x y z coords + geometry = conformer + calculation = 'opt' + torsion_data = f' \n' + + if Descriptors.NumRadicalElectrons(mol) == 0: + mult_chg = '0 1' + if Descriptors.NumRadicalElectrons(mol) == 1: + mult_chg = '0 2' + if Descriptors.NumRadicalElectrons(mol) == 2: + mult_chg = '0 3' + + file_name = f'{mol_name}_{job_name}_xTB.py' + #Includes information about basis set to allow for ramping + title = f'{mol_name} {job_name} Smile String: {smile}' + number_proc = 'set_num_threads(4)' + memory_num = '4000 mb' + number_of_atoms = mol.GetNumAtoms() + + with open(file_name, 'w') as file: + file.write('import subprocess\n') + file.write('import tempfile\n') + file.write('import os\n') + file.write('import psi4\n') + file.write('from QCflow.write_xTB import *\n') + file.write(' \n')# + file.write(f'xyz_string = """{number_of_atoms}\n') + file.write(' \n') + if old_chk==' ': + for atom,symbol in enumerate(symbols): + p = geometry.GetAtomPosition(atom) + # atom x y z + line = f' {symbol} {p.x:.5f} {p.y:.5f} {p.z:.5f} \n' + file.write(line) + file.write('""" \n') + file.write(' \n') + file.write('run_xtb_optimisation(xyz_string, "optimised_structure.xyz", gfn=2)\n') + file.write(' \n') + file.write(f"mol_sp = psi4.geometry(open('optimised_structure.xyz').read())\n")#use geometry from opt + file.write('psi4.core.set_active_molecule(mol_sp)\n')#set as active molecule + file.write(f"mol_sp.set_molecular_charge(0) \n")# set charge + file.write(f"mol_sp.set_multiplicity(1) \n")# set multiplicity + file.write(' \n') + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df'}})\n") + file.write(' \n') + file.write(f"energy, wfn = psi4.energy('{functional}', return_wfn=True) \n") + file.write(' \n') + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('hartree_to_ev = 27.2114079527\n') + file.write('energy_ev = energy * hartree_to_ev\n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_xTB_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Single Point energy: {energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + file.write(' \n') + + +def write_xTB_psi4_reorg(job_name, mol_name, functional='b3lyp', basis_set='6-31g*'): + """ + Generates a xTB and psi4 input file based on the provided parameters for reorganisation calculations. + + Parameters + ---------- + job_name (str): The type of job to run. Possible values: + - 'cation': Geometry optimisation cation (opt_c) and single of neutral charge, cation geometry (n_c_geo) + - 'anion': Geometry optimisation anion (opt_a) and single of neutral charge, anion geometry (n_a_geo) + - 'sp_c': Single point calculation of neutral geometry at cation charge + - 'sp_a': Single point calculation of neutral geometry at anion charge + mol_name (str): The name of the molecule. + functional (str, optional): The functional to use. Default is 'b3lyp'. + basis_set (str, optional): The basis set to use. Default is '6-31g'. + + Returns + ------- + None: Writes the psi4 input file to disk. + + """ + if (job_name=='cation'): + + charge = 1 + mult = 2 + + job_type = 'opt_c' + + job_sp = 'n_c_geo' + + with open('optimised_structure.xyz', 'r') as f: + lines = f.readlines() + number_of_atoms = int(lines[0].strip()) + coordinate_lines = [line.strip() for line in lines[2:2 + number_of_atoms]] + + + if (job_name=='anion'): + + charge = -1 + mult = 2 + + job_type = 'opt_a' + + job_sp = 'n_a_geo' + + with open('optimised_structure.xyz', 'r') as f: + lines = f.readlines() + number_of_atoms = int(lines[0].strip()) + coordinate_lines = [line.strip() for line in lines[2:2 + number_of_atoms]] + + if (job_name=='sp_c'): + + charge = 1 + mult = 2 + + with open('optimised_structure.xyz', 'r') as f: + lines = f.readlines() + number_of_atoms = int(lines[0].strip()) + coordinate_lines = [line.strip() for line in lines[2:2 + number_of_atoms]] + + if (job_name=='sp_a'): + + charge = -1 + mult = 2 + + with open('optimised_structure.xyz', 'r') as f: + lines = f.readlines() + number_of_atoms = int(lines[0].strip()) + coordinate_lines = [line.strip() for line in lines[2:2 + number_of_atoms]] + + + file_name = f'{mol_name}_{job_name}.py' + #Includes information about basis set to allow for ramping + memory_num = '4000 mb' + + if (job_name=='anion') or (job_name=='cation'): + + with open(file_name, 'w') as file: + file.write('import subprocess\n') + file.write('import tempfile\n') + file.write('import os\n') + file.write('import psi4\n') + file.write('from QCflow.write_xTB import *\n') + file.write(' \n')# + file.write(f'xyz_string = """{number_of_atoms}\n') + file.write(' \n') + # Now insert the parsed coordinate lines + for line in coordinate_lines: + file.write(f'{line}\n') + file.write('""" \n') + file.write(' \n') + file.write(f"run_xtb_optimisation(xyz_string, '{job_type}_optimised_structure.xyz', charge={charge}, uhf={mult}, gfn=2)\n") + file.write(' \n') + file.write(f"mol = psi4.geometry(open('{job_type}_optimised_structure.xyz').read())\n")#use geometry from opt + file.write('psi4.core.set_active_molecule(mol)\n')#set as active molecule + file.write(f"mol.set_molecular_charge({charge}) \n")# set charge + file.write(f"mol.set_multiplicity({mult}) \n")# set multiplicity + file.write("psi4.set_module_options('scf', {'maxiter': 100}) \n") + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df', 'reference': 'uhf'}})\n") + file.write(' \n') + file.write(f"energy, wfn = psi4.optimize('{functional}', engine='geometric', return_wfn=True) \n") + file.write(' \n') + file.write("optimized_geometry_xyz = wfn.molecule().to_string('xyz') \n") + file.write(f"with open('{mol_name}_{job_type}.xyz', 'w') as xyz_file:\n") + file.write(" xyz_file.write(optimized_geometry_xyz) \n") + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('hartree_to_ev = 27.2114079527 \n') + file.write('energy_ev = energy * hartree_to_ev\n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_{job_type}_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Optimized energy: {energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + file.write(' \n')# Now do the n_c_geo or n_a_geo part + file.write(f"mol_sp = psi4.geometry(open('{mol_name}_{job_type}.xyz').read())\n")#use geometry from opt + file.write('psi4.core.set_active_molecule(mol_sp)\n')#set as active molecule + file.write(f"mol_sp.set_molecular_charge(0) \n")# set charge + file.write(f"mol_sp.set_multiplicity(1) \n")# set multiplicity + file.write(' \n') + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df'}})\n") + file.write(' \n') + file.write(f"energy, wfn = psi4.energy('{functional}', return_wfn=True) \n") + file.write(' \n') + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('energy_ev = energy * hartree_to_ev\n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_{job_sp}_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Single Point energy: {energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + file.write(' \n') + + if (job_name=='sp_c') or (job_name=='sp_a'): + with open(file_name, 'w') as file: + file.write('import psi4\n') + file.write(' \n')# + file.write(f"psi4.set_memory('{memory_num}')\n") + file.write('psi4.set_num_threads(4)\n') + file.write(' \n')# + file.write(f"mol = psi4.geometry(open('optimised_structure.xyz').read())\n")#use geometry from opt + file.write('psi4.core.set_active_molecule(mol)\n')#set as active molecule + file.write(f"mol.set_molecular_charge({charge}) \n")# set charge + file.write(f"mol.set_multiplicity({mult}) \n")# set multiplicity + file.write(' \n') + file.write("psi4.set_module_options('scf', {'maxiter': 50}) \n") + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df', 'reference': 'uhf'}})\n") + file.write(' \n') + file.write(f"energy, wfn = psi4.energy('{functional}', return_wfn=True) \n") + file.write(' \n') + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('hartree_to_ev = 27.2114079527 \n') + file.write('energy_ev = energy * hartree_to_ev\n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_{job_name}_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Single Point energy: {energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + file.write(' \n') diff --git a/README.md b/README.md index 2d76b34..87665e2 100755 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![QCflow stable](https://github.com/matta-research-group/QCflow/actions/workflows/run_test.yml/badge.svg?branch=qcflow-psi4)](https://github.com/matta-research-group/QCflow/actions/workflows/run_test.yml) ![Python](https://img.shields.io/badge/language-Python-blue.svg) [![Documentation](https://img.shields.io/badge/Documentation-Online-brightgreen)](https://matta-research-group.github.io/QCflow/) -[![GitHub Last commit](https://img.shields.io/github/last-commit/matta-research-group/QCflow)](https://github.com/matta-research-group/QCflow/commits/qcflow-0.3/) +[![GitHub Last commit](https://img.shields.io/github/last-commit/matta-research-group/QCflow)](https://github.com/matta-research-group/QCflow/commits/qcflow-0.5/) [![GitHub stars](https://img.shields.io/github/stars/matta-research-group/QCflow)](https://github.com/matta-research-group/QCflow/stargazers) A cheminformatics -> quantum chemistry workflow toolkit leveraging [rdkit](https://github.com/rdkit/rdkit) and [cclib](https://github.com/cclib/cclib). @@ -120,6 +120,7 @@ QCflow can prepare and submit input files for the following Gaussian 16 jobs: QCflow can prepare and submit input files for the following Psi4 jobs: - Single point calculation, neutral → `sp` - Geometry optimisation, neutral → `opt` +- Geometry optimisation, neutral from txt file geometry → `opt_pre_geom` - Geometry optimisation cation (opt_c) and single of neutral charge, cation geometry (n_c_geo) → `cation` - Geometry optimisation anion (opt_a) and single of neutral charge, anion geometry (n_a_geo) → `anion` - Single point calculation, anion → `sp_a` @@ -135,80 +136,9 @@ Additionally, we welcome any developments and improvements user may have and if ```bash . -├── Examples -│   ├── example_calculation_submission.py -│   ├── example_data_extraction.ipynb -│   ├── example_data_extraction_bulk.py -│   ├── example_molecule_building.ipynb -│   └── example_molecule_logs -│   ├── 18 -│   │   └── 18_opt.log -│   ├── b -│   │   └── b_opt.log -│   └── b_18_single_v2 -│   ├── b_18_single_v2_n_a_geo.log -│   ├── b_18_single_v2_n_c_geo.log -│   ├── b_18_single_v2_opt.log -│   ├── b_18_single_v2_opt_a.log -│   ├── b_18_single_v2_opt_c.log -│   ├── b_18_single_v2_sp_a.log -│   └── b_18_single_v2_sp_c.log -├── LICENSE -├── QCflow -│   ├── __init__.py -│   ├── __pycache__ -│   │   ├── __init__.cpython-311.pyc -│   │   ├── __init__.cpython-313.pyc -│   │   ├── __init__.cpython-38.pyc -│   │   ├── __init__.cpython-39.pyc -│   │   ├── energy_calculations.cpython-313.pyc -│   │   ├── find_torsion.cpython-311.pyc -│   │   ├── find_torsion.cpython-313.pyc -│   │   ├── find_torsion.cpython-39.pyc -│   │   ├── fragments.cpython-311.pyc -│   │   ├── fragments.cpython-313.pyc -│   │   ├── fragments.cpython-38.pyc -│   │   ├── fragments.cpython-39.pyc -│   │   ├── load_gaussian.cpython-311.pyc -│   │   ├── load_gaussian.cpython-313.pyc -│   │   ├── load_gaussian.cpython-39.pyc -│   │   ├── run_gaussian.cpython-313.pyc -│   │   ├── run_opt_neutral.cpython-311.pyc -│   │   ├── run_opt_neutral.cpython-39.pyc -│   │   ├── run_opt_set_dihedral.cpython-311.pyc -│   │   ├── run_psi4.cpython-313.pyc -│   │   ├── sa_score.cpython-313.pyc -│   │   ├── slurm.cpython-311.pyc -│   │   ├── slurm.cpython-313.pyc -│   │   ├── slurm.cpython-39.pyc -│   │   ├── testing_data.cpython-313.pyc -│   │   ├── torsion_parser.cpython-311.pyc -│   │   ├── torsion_parser.cpython-313.pyc -│   │   ├── torsion_parser.cpython-39.pyc -│   │   ├── torsion_run.cpython-311.pyc -│   │   ├── torsion_run.cpython-39.pyc -│   │   ├── write_gaussian.cpython-311.pyc -│   │   ├── write_gaussian.cpython-313.pyc -│   │   ├── write_gaussian.cpython-39.pyc -│   │   └── write_psi4.cpython-313.pyc -│   ├── energy_calculations.py -│   ├── find_torsion.py -│   ├── fragments.py -│   ├── future_functions -│   │   └── orbital_parse.py -│   ├── load_gaussian.py -│   ├── run_gaussian.py -│   ├── run_psi4.py -│   ├── sa_score.py -│   ├── slurm.py -│   ├── testing_data.py -│   ├── torsion_parser.py -│   ├── write_gaussian.py -│   └── write_psi4.py -├── QCflow.yml -├── QCflow_logo_narrow.jpg -├── README.md ├── build +│   ├── bdist.linux-x86_64 +│   ├── bdist.macosx-11.0-arm64 │   └── lib │   └── QCflow │   ├── __init__.py @@ -218,15 +148,15 @@ Additionally, we welcome any developments and improvements user may have and if │   ├── load_gaussian.py │   ├── run_gaussian.py │   ├── run_psi4.py +│   ├── run_xTB.py │   ├── sa_score.py │   ├── slurm.py │   ├── testing_data.py │   ├── torsion_parser.py │   ├── write_gaussian.py -│   └── write_psi4.py +│   ├── write_psi4.py +│   └── write_xTB.py ├── docs -│   ├── Makefile -│   ├── README.md │   ├── _build │   │   ├── doctrees │   │   │   ├── api.doctree @@ -235,35 +165,30 @@ Additionally, we welcome any developments and improvements user may have and if │   │   │   └── index.doctree │   │   └── html │   │   ├── _modules -│   │   │   ├── QCflow -│   │   │   │   ├── energy_calculations.html -│   │   │   │   ├── find_torsion.html -│   │   │   │   ├── fragments.html -│   │   │   │   ├── load_gaussian.html -│   │   │   │   ├── run_gaussian.html -│   │   │   │   ├── run_psi4.html -│   │   │   │   ├── sa_score.html -│   │   │   │   ├── slurm.html -│   │   │   │   ├── testing_data.html -│   │   │   │   ├── torsion_parser.html -│   │   │   │   ├── write_gaussian.html -│   │   │   │   └── write_psi4.html -│   │   │   └── index.html +│   │   │   ├── index.html +│   │   │   └── QCflow +│   │   │   ├── energy_calculations.html +│   │   │   ├── find_torsion.html +│   │   │   ├── fragments.html +│   │   │   ├── load_gaussian.html +│   │   │   ├── run_gaussian.html +│   │   │   ├── run_psi4.html +│   │   │   ├── sa_score.html +│   │   │   ├── slurm.html +│   │   │   ├── testing_data.html +│   │   │   ├── torsion_parser.html +│   │   │   ├── write_gaussian.html +│   │   │   └── write_psi4.html │   │   ├── _sources │   │   │   ├── api.rst.txt │   │   │   ├── getting_started.rst.txt │   │   │   └── index.rst.txt │   │   ├── _static -│   │   │   ├── README.md │   │   │   ├── _sphinx_javascript_frameworks_compat.js │   │   │   ├── basic.css │   │   │   ├── css │   │   │   │   ├── badge_only.css │   │   │   │   ├── fonts -│   │   │   │   │   ├── Roboto-Slab-Bold.woff -│   │   │   │   │   ├── Roboto-Slab-Bold.woff2 -│   │   │   │   │   ├── Roboto-Slab-Regular.woff -│   │   │   │   │   ├── Roboto-Slab-Regular.woff2 │   │   │   │   │   ├── fontawesome-webfont.eot │   │   │   │   │   ├── fontawesome-webfont.svg │   │   │   │   │   ├── fontawesome-webfont.ttf @@ -276,7 +201,11 @@ Additionally, we welcome any developments and improvements user may have and if │   │   │   │   │   ├── lato-normal-italic.woff │   │   │   │   │   ├── lato-normal-italic.woff2 │   │   │   │   │   ├── lato-normal.woff -│   │   │   │   │   └── lato-normal.woff2 +│   │   │   │   │   ├── lato-normal.woff2 +│   │   │   │   │   ├── Roboto-Slab-Bold.woff +│   │   │   │   │   ├── Roboto-Slab-Bold.woff2 +│   │   │   │   │   ├── Roboto-Slab-Regular.woff +│   │   │   │   │   └── Roboto-Slab-Regular.woff2 │   │   │   │   └── theme.css │   │   │   ├── doctools.js │   │   │   ├── documentation_options.js @@ -317,6 +246,7 @@ Additionally, we welcome any developments and improvements user may have and if │   │   │   ├── minus.png │   │   │   ├── plus.png │   │   │   ├── pygments.css +│   │   │   ├── README.md │   │   │   ├── searchtools.js │   │   │   └── sphinx_highlight.js │   │   ├── api.html @@ -336,15 +266,102 @@ Additionally, we welcome any developments and improvements user may have and if │   ├── getting_started.rst │   ├── index.rst │   ├── make.bat +│   ├── Makefile +│   ├── README.md │   ├── requirements.yaml │   └── timer.dat +├── Examples +│   ├── example_data_extraction_bulk.py +│   ├── example_data_extraction.ipynb +│   ├── example_g16_calculation_submission.py +│   ├── example_molecule_building.ipynb +│   ├── example_molecule_logs +│   │   ├── 18 +│   │   │   └── 18_opt.log +│   │   ├── b +│   │   │   └── b_opt.log +│   │   └── b_18_single_v2 +│   │   ├── b_18_single_v2_n_a_geo.log +│   │   ├── b_18_single_v2_n_c_geo.log +│   │   ├── b_18_single_v2_opt_a.log +│   │   ├── b_18_single_v2_opt_c.log +│   │   ├── b_18_single_v2_opt.log +│   │   ├── b_18_single_v2_sp_a.log +│   │   └── b_18_single_v2_sp_c.log +│   ├── example_psi4_opt_submission.py +│   ├── example_torsional_scan_script.py +│   └── FilteredSmilesFragments_props.csv +├── LICENSE +├── QCflow +│   ├── __init__.py +│   ├── __pycache__ +│   │   ├── __init__.cpython-311.pyc +│   │   ├── __init__.cpython-313.pyc +│   │   ├── __init__.cpython-38.pyc +│   │   ├── __init__.cpython-39.pyc +│   │   ├── energy_calculations.cpython-313.pyc +│   │   ├── find_torsion.cpython-311.pyc +│   │   ├── find_torsion.cpython-313.pyc +│   │   ├── find_torsion.cpython-39.pyc +│   │   ├── fragments.cpython-311.pyc +│   │   ├── fragments.cpython-313.pyc +│   │   ├── fragments.cpython-38.pyc +│   │   ├── fragments.cpython-39.pyc +│   │   ├── load_gaussian.cpython-311.pyc +│   │   ├── load_gaussian.cpython-313.pyc +│   │   ├── load_gaussian.cpython-39.pyc +│   │   ├── run_gaussian.cpython-313.pyc +│   │   ├── run_opt_neutral.cpython-311.pyc +│   │   ├── run_opt_neutral.cpython-39.pyc +│   │   ├── run_opt_set_dihedral.cpython-311.pyc +│   │   ├── run_psi4.cpython-313.pyc +│   │   ├── run_xTB.cpython-313.pyc +│   │   ├── sa_score.cpython-313.pyc +│   │   ├── slurm.cpython-311.pyc +│   │   ├── slurm.cpython-313.pyc +│   │   ├── slurm.cpython-39.pyc +│   │   ├── testing_data.cpython-313.pyc +│   │   ├── torsion_parser.cpython-311.pyc +│   │   ├── torsion_parser.cpython-313.pyc +│   │   ├── torsion_parser.cpython-39.pyc +│   │   ├── torsion_run.cpython-311.pyc +│   │   ├── torsion_run.cpython-39.pyc +│   │   ├── write_gaussian.cpython-311.pyc +│   │   ├── write_gaussian.cpython-313.pyc +│   │   ├── write_gaussian.cpython-39.pyc +│   │   ├── write_psi4.cpython-313.pyc +│   │   └── write_xTB.cpython-313.pyc +│   ├── energy_calculations.py +│   ├── find_torsion.py +│   ├── fragments.py +│   ├── future_functions +│   │   └── orbital_parse.py +│   ├── load_gaussian.py +│   ├── run_gaussian.py +│   ├── run_psi4.py +│   ├── run_xTB.py +│   ├── sa_score.py +│   ├── slurm.py +│   ├── testing_data.py +│   ├── torsion_parser.py +│   ├── write_gaussian.py +│   ├── write_psi4.py +│   └── write_xTB.py +├── QCflow_logo_narrow.jpg ├── qcflow.egg-info +│   ├── dependency_links.txt │   ├── PKG-INFO │   ├── SOURCES.txt -│   ├── dependency_links.txt │   └── top_level.txt +├── QCflow.yml +├── README.md ├── setup.py └── tests + ├── __pycache__ + │   ├── test.cpython-313-pytest-7.4.4.pyc + │   ├── test.cpython-313-pytest-8.3.3.pyc + │   ├── tests.cpython-312-pytest-7.4.4.pyc + │   └── tests.cpython-312.pyc ├── 1 │   ├── 1_cation.py │   ├── 1_opt.py @@ -360,26 +377,26 @@ Additionally, we welcome any developments and improvements user may have and if │   ├── 12_cation.sh │   ├── 12_n_a_geo_energy_and_gap.txt │   ├── 12_n_c_geo_energy_and_gap.txt + │   ├── 12_opt_a_energy_and_gap.txt + │   ├── 12_opt_a.xyz + │   ├── 12_opt_c_energy_and_gap.txt + │   ├── 12_opt_c.xyz + │   ├── 12_opt_energy_and_gap.txt │   ├── 12_opt.err │   ├── 12_opt.out │   ├── 12_opt.py │   ├── 12_opt.sh │   ├── 12_opt.xyz - │   ├── 12_opt_a.xyz - │   ├── 12_opt_a_energy_and_gap.txt - │   ├── 12_opt_c.xyz - │   ├── 12_opt_c_energy_and_gap.txt - │   ├── 12_opt_energy_and_gap.txt + │   ├── 12_sp_a_energy_and_gap.txt │   ├── 12_sp_a.err │   ├── 12_sp_a.out │   ├── 12_sp_a.py │   ├── 12_sp_a.sh - │   ├── 12_sp_a_energy_and_gap.txt + │   ├── 12_sp_c_energy_and_gap.txt │   ├── 12_sp_c.err │   ├── 12_sp_c.out │   ├── 12_sp_c.py │   ├── 12_sp_c.sh - │   ├── 12_sp_c_energy_and_gap.txt │   └── timer.dat ├── 145 │   ├── 145_anion.err @@ -392,48 +409,43 @@ Additionally, we welcome any developments and improvements user may have and if │   ├── 145_cation.sh │   ├── 145_n_a_geo_energy_and_gap.txt │   ├── 145_n_c_geo_energy_and_gap.txt + │   ├── 145_opt_a_energy_and_gap.txt + │   ├── 145_opt_a.xyz + │   ├── 145_opt_c_energy_and_gap.txt + │   ├── 145_opt_c.xyz + │   ├── 145_opt_energy_and_gap.txt │   ├── 145_opt.err │   ├── 145_opt.out │   ├── 145_opt.py │   ├── 145_opt.sh │   ├── 145_opt.xyz - │   ├── 145_opt_a.xyz - │   ├── 145_opt_a_energy_and_gap.txt - │   ├── 145_opt_c.xyz - │   ├── 145_opt_c_energy_and_gap.txt - │   ├── 145_opt_energy_and_gap.txt + │   ├── 145_sp_a_energy_and_gap.txt │   ├── 145_sp_a.err │   ├── 145_sp_a.out │   ├── 145_sp_a.py │   ├── 145_sp_a.sh - │   ├── 145_sp_a_energy_and_gap.txt + │   ├── 145_sp_c_energy_and_gap.txt │   ├── 145_sp_c.err │   ├── 145_sp_c.out │   ├── 145_sp_c.py │   ├── 145_sp_c.sh - │   ├── 145_sp_c_energy_and_gap.txt │   └── timer.dat - ├── __pycache__ - │   ├── test.cpython-313-pytest-7.4.4.pyc - │   ├── test.cpython-313-pytest-8.3.3.pyc - │   ├── tests.cpython-312-pytest-7.4.4.pyc - │   └── tests.cpython-312.pyc ├── b_18_v2 │   ├── b_18_v2_n_a_geo.log │   ├── b_18_v2_n_c_geo.log + │   ├── b_18_v2_opt_a.log + │   ├── b_18_v2_opt_c.log │   ├── b_18_v2_opt.com │   ├── b_18_v2_opt.log │   ├── b_18_v2_opt.sh - │   ├── b_18_v2_opt_a.log - │   ├── b_18_v2_opt_c.log │   ├── b_18_v2_sp_a.log │   ├── b_18_v2_sp_c.log │   ├── b_18_v2_sp_hirsh.log │   └── b_18_v2_tor.log ├── example_dic.json ├── fake_dict.json - ├── test.py ├── test_dict.json + ├── test.py ├── timer.dat └── torsion_1 └── torsion_1_tor.log diff --git a/build/lib/QCflow/__init__.py b/build/lib/QCflow/__init__.py index 72ea573..77dbda4 100644 --- a/build/lib/QCflow/__init__.py +++ b/build/lib/QCflow/__init__.py @@ -10,3 +10,5 @@ from . import write_psi4 from . import run_psi4 from . import sa_score +from . import write_xTB +from . import run_xTB diff --git a/build/lib/QCflow/run_psi4.py b/build/lib/QCflow/run_psi4.py index 417138f..7009f82 100644 --- a/build/lib/QCflow/run_psi4.py +++ b/build/lib/QCflow/run_psi4.py @@ -17,6 +17,7 @@ def run_psi4(job_name, mol_name, mol_smile, time=24, cpus=10, functional='b3lyp' job_name (str): The type of job to run. Possible values: - 'sp': Single Point neutral - 'opt': Optimisation neutral + - 'opt_pre_geom': Optimisation neutral where a txt file called '{mol_name}_opt.xyz' is present in the molecule directory. This file should contain the geometry to be used for the optimisation in XYZ format. If this file is not present, the function will default to 'opt'. - 'cation': Geometry optimisation cation (opt_c) and single of neutral charge, cation geometry (n_c_geo) - 'anion': Geometry optimisation anion (opt_a) and single of neutral charge, anion geometry (n_a_geo) - 'sp_c': Single point calculation of neutral geometry at cation charge @@ -50,7 +51,14 @@ def run_psi4(job_name, mol_name, mol_smile, time=24, cpus=10, functional='b3lyp' os.mkdir(f'{mol_name}') #makes a directory for the molecule os.chdir(f'{mol_name}') #goes into directory - if (job_name=='opt') or (job_name=='sp'): + if (job_name=='opt') or (job_name=='sp') or (job_name=='opt_pre_geom'): + + if (job_name=='opt_pre_geom'): + if os.path.exists(f'{mol_name}_opt.xyz'): + print(f"Found {mol_name}_opt.xyz, using this geometry for optimisation.") + else: + print(f"{mol_name}_opt.xyz not found, defaulting to 'opt' job type.") + job_name = 'opt' #turns smiles string into rdkit object mol = Chem.MolFromSmiles(mol_smile) #gets rdkit estimated coordinates of dimer @@ -64,6 +72,9 @@ def run_psi4(job_name, mol_name, mol_smile, time=24, cpus=10, functional='b3lyp' if (job_name=='cation') or (job_name=='anion') or (job_name=='sp_c') or (job_name=='sp_a'): write_psi4_reorg(job_name, mol_name, functional, basis_set) + + else: + print(f"Invalid job_name: {job_name}, please use 'opt', 'sp', 'opt_pre_geom', 'cation', 'anion', 'sp_c' or 'sp_a'. See documentation for more details.") #writes the slurm file diff --git a/build/lib/QCflow/run_xTB.py b/build/lib/QCflow/run_xTB.py new file mode 100644 index 0000000..0a88d28 --- /dev/null +++ b/build/lib/QCflow/run_xTB.py @@ -0,0 +1,76 @@ +from QCflow.fragments import * +from QCflow.find_torsion import * +from QCflow.write_gaussian import * +from QCflow.slurm import * +from QCflow.torsion_parser import * +from QCflow.run_gaussian import * +from QCflow.write_psi4 import * +from QCflow.write_xTB import * +import os +import json +import subprocess +import tempfile + +def run_xTB(job_name, mol_name, mol_smile, time=2, cpus=10, functional='b3lyp', basis_set='6-31g*'): + """ + Submits a xTB calculation to CREATE HPC. + + Parameters + ---------- + job_name (str): The type of job to run. Possible values: + - 'opt': Optimisation neutral + - 'cation': Geometry optimisation cation (opt_c) and single of neutral charge, cation geometry (n_c_geo) + - 'anion': Geometry optimisation anion (opt_a) and single of neutral charge, anion geometry (n_a_geo) + - 'sp_c': Single point calculation of neutral geometry at cation charge + - 'sp_a': Single point calculation of neutral geometry at anion charge + mol_name : str + Name of the oligomer. + mol_smile : str + SMILES string of the oligomer. + time (int, optional): The time limit for the job in hours. Default is 24. (Max is 48) + cpus (int, optional): The number of CPUs to allocate for the job. Default is 10. + functional : str, optional + Quantum chemistry functional to be used (default is 'B3LYP'). + basis_set : str, optional + Basis set to be used (default is '6-31G*'). + + + Returns + ------- + - For 'opt': + - Converts the SMILES string to an RDKit molecule object. + - Generates 3D coordinates for the molecule. + - Predicts the conformer geometry. + - Runs xTB optimisation using the generated coordinates. + - Writes a psi4 input file with the optimised geometry. + - Runs a single point calculation. + Finally, the function writes a SLURM script, submits the job, and returns to the previous directory. + """ + if os.path.exists(f'{mol_name}'): + os.chdir(f'{mol_name}') #goes into directory + else: + os.mkdir(f'{mol_name}') #makes a directory for the molecule + os.chdir(f'{mol_name}') #goes into directory + + if (job_name=='opt'): + #turns smiles string into rdkit object + mol = Chem.MolFromSmiles(mol_smile) + #gets rdkit estimated coordinates of dimer + mol3d = embed_molecule(mol) + + conf_geo = rdkit_predict_conf(mol_smile) + #writes a guassian input file + write_xTB_psi4(job_name, mol_name, mol_smile, functional, basis_set, mol=mol3d, conformer=conf_geo) + + #for reorganisation calcultions + if (job_name=='cation') or (job_name=='anion') or (job_name=='sp_c') or (job_name=='sp_a'): + + write_xTB_psi4_reorg(job_name, mol_name, functional, basis_set) + + + #writes the slurm file + write_slurm_psi4(job_name, mol_name, time, cpus) + #submits the slurm jon + submit_slurm_job(job_name, mol_name) + #goes back to previous dirctory + os.chdir(os.path.dirname(os.getcwd())) diff --git a/build/lib/QCflow/write_psi4.py b/build/lib/QCflow/write_psi4.py old mode 100644 new mode 100755 index 322f9b0..92c8746 --- a/build/lib/QCflow/write_psi4.py +++ b/build/lib/QCflow/write_psi4.py @@ -20,7 +20,8 @@ def write_psi4(job_name, mol_name, smile, functional='b3lyp', basis_set='6-31g*' ---------- job_name (str): The type of job to run. Possible values: - 'sp': Single Point neutral - - 'opt': Optimisation neutral + - 'opt': Optimisation neutral + - 'opt_pre_geom' : Optimisation neutral, from stated geometry, where a txt file called '{mol_name}_opt.xyz' is present and in psi4 format mol_name (str): The name of the molecule. smile (str): The SMILE string of the molecule. functional (str, optional): The functional to use. Default is 'b3lyp'. @@ -66,7 +67,29 @@ def write_psi4(job_name, mol_name, smile, functional='b3lyp', basis_set='6-31g*' if Descriptors.NumRadicalElectrons(mol) == 1: mult_chg = '0 2' if Descriptors.NumRadicalElectrons(mol) == 2: - mult_chg = '0 3' + mult_chg = '0 3' + + if (job_name=='opt_pre_geom'): + + old_chk = ' ' + + geometry_file = f'{mol_name}_opt.xyz' + calculation = 'opt' + torsion_data = f' \n' + + if Descriptors.NumRadicalElectrons(mol) == 0: + charge = '0' + mult = '1' + mult_chg = f'{charge} {mult}' + if Descriptors.NumRadicalElectrons(mol) == 1: + charge = '0' + mult = '2' + mult_chg = f'{charge} {mult}' + if Descriptors.NumRadicalElectrons(mol) == 2: + charge = '0' + mult = '3' + mult_chg = f'{charge} {mult}' + file_name = f'{mol_name}_{job_name}.py' #Includes information about basis set to allow for ramping @@ -74,56 +97,98 @@ def write_psi4(job_name, mol_name, smile, functional='b3lyp', basis_set='6-31g*' number_proc = 'set_num_threads(4)' memory_num = '4000 mb' - - with open(file_name, 'w') as file: - file.write('import psi4\n') - file.write(' \n')# - file.write(f"psi4.set_memory('{memory_num}')\n") - file.write('psi4.set_num_threads(4)\n') - file.write(' \n')# - file.write('mol_name = psi4.geometry("""\n')# - #file.write(' symmetry c1 \n') - file.write(f' {mult_chg} \n')# - # if no previous checkpoint file, get geometry - if old_chk==' ': - for atom,symbol in enumerate(symbols): - p = geometry.GetAtomPosition(atom) - # atom x y z - line = f' {symbol} {p.x:.5f} {p.y:.5f} {p.z:.5f} \n' - file.write(line) - file.write('""") \n') - file.write(' \n') - file.write("psi4.set_module_options('scf', {'maxiter': 75}) \n") - file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df'}})\n") - file.write(' \n') - if (job_name=='sp'): - file.write(f"energy, wfn = psi4.energy('{functional}', return_wfn=True) \n") - file.write(' \n') - file.write("sp_geometry_xyz = wfn.molecule().to_string('xyz') \n") - file.write(f"with open('{mol_name}_{job_name}.xyz', 'w') as xyz_file:\n") - file.write(" xyz_file.write(sp_geometry_xyz) \n") - file.write(' \n') - file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") - file.write(' \n') - file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') - file.write(' \n') - file.write('n_occ_alpha = wfn.nalpha() \n') - file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') - file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') - file.write('energy_gap = lumo_energy - homo_energy \n') + if (job_name== 'sp') or (job_name=='opt'): + with open(file_name, 'w') as file: + file.write('import psi4\n') + file.write(' \n')# + file.write(f"psi4.set_memory('{memory_num}')\n") + file.write('psi4.set_num_threads(4)\n') + file.write(' \n')# + file.write('mol_name = psi4.geometry("""\n')# + #file.write(' symmetry c1 \n') + file.write(f' {mult_chg} \n')# + # if no previous checkpoint file, get geometry + if old_chk==' ': + for atom,symbol in enumerate(symbols): + p = geometry.GetAtomPosition(atom) + # atom x y z + line = f' {symbol} {p.x:.5f} {p.y:.5f} {p.z:.5f} \n' + file.write(line) + file.write('""") \n') + file.write(' \n') + file.write("psi4.set_options({'SCF__MAXITER': 100, 'GEOM_MAXITER': 200}) \n") + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df'}})\n") file.write(' \n') - file.write('hartree_to_ev = 27.2114079527 \n') - file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') - file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') - file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + if (job_name=='sp'): + file.write(f"energy, wfn = psi4.energy('{functional}', return_wfn=True) \n") + file.write(' \n') + file.write("sp_geometry_xyz = wfn.molecule().to_string('xyz') \n") + file.write(f"with open('{mol_name}_{job_name}.xyz', 'w') as xyz_file:\n") + file.write(" xyz_file.write(sp_geometry_xyz) \n") + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('hartree_to_ev = 27.2114079527 \n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_{job_name}_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Optimized energy: {energy:.6f} Hatree\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + if (job_name=='opt'): + file.write(f"energy, wfn = psi4.optimize('{functional}', engine='geometric', return_wfn=True) \n") + file.write(' \n') + file.write("optimized_geometry_xyz = wfn.molecule().to_string('xyz') \n") + file.write(f"with open('{mol_name}_{job_name}.xyz', 'w') as xyz_file:\n") + file.write(" xyz_file.write(optimized_geometry_xyz) \n") + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('hartree_to_ev = 27.2114079527 \n') + file.write('energy_ev = energy * hartree_to_ev \n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_{job_name}_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Optimized energy: {energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + file.write(' \n') + + if (job_name == 'opt_pre_geom'): + with open(file_name, 'w') as file: + file.write('import psi4\n') + file.write(' \n')# + file.write(f"psi4.set_memory('{memory_num}')\n") + file.write('psi4.set_num_threads(4)\n') + file.write(' \n')# + file.write(f"mol = psi4.geometry(open('{geometry_file}').read())\n")#use geometry from opt + file.write('psi4.core.set_active_molecule(mol)\n')#set as active molecule + file.write(f"mol.set_molecular_charge({charge}) \n")# set charge + file.write(f"mol.set_multiplicity({mult}) \n")# set multiplicity + file.write("psi4.set_options({'SCF__MAXITER': 100, 'GEOM_MAXITER': 200}) \n") + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df'}})\n") file.write(' \n') - file.write(f"with open('{mol_name}_{job_name}_energy_and_gap.txt', 'w') as file:\n") - file.write(' file.write(f"Optimized energy: {energy:.6f} Hatree\\n") \n') - file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') - file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') - file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') - if (job_name=='opt'): - file.write(f"energy, wfn = psi4.optimize('{functional}', return_wfn=True) \n") + file.write(f"energy, wfn = psi4.optimize('{functional}', engine='geometric', return_wfn=True) \n") file.write(' \n') file.write("optimized_geometry_xyz = wfn.molecule().to_string('xyz') \n") file.write(f"with open('{mol_name}_{job_name}.xyz', 'w') as xyz_file:\n") @@ -149,7 +214,10 @@ def write_psi4(job_name, mol_name, smile, functional='b3lyp', basis_set='6-31g*' file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') - file.write(' \n') + file.write(' \n') + + else: + print(f"Invalid job_name: {job_name}, please use 'opt', 'sp' or 'opt_pre_geom'. See documentation for more details.") def write_psi4_reorg(job_name, mol_name, functional='b3lyp', basis_set='6-31g*'): @@ -227,10 +295,10 @@ def write_psi4_reorg(job_name, mol_name, functional='b3lyp', basis_set='6-31g*') file.write(f"mol.set_multiplicity({mult}) \n")# set multiplicity file.write(' \n') if (job_name=='cation') or (job_name=='anion'): - file.write("psi4.set_module_options('scf', {'maxiter': 100}) \n") + file.write("psi4.set_options({'SCF__MAXITER': 100, 'GEOM_MAXITER': 200}) \n") file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df', 'reference': 'uhf'}})\n") file.write(' \n') - file.write(f"energy, wfn = psi4.optimize('{functional}', return_wfn=True) \n") + file.write(f"energy, wfn = psi4.optimize('{functional}', engine='geometric', return_wfn=True) \n") file.write(' \n') file.write("optimized_geometry_xyz = wfn.molecule().to_string('xyz') \n") file.write(f"with open('{mol_name}_{job_type}.xyz', 'w') as xyz_file:\n") @@ -316,5 +384,7 @@ def write_psi4_reorg(job_name, mol_name, functional='b3lyp', basis_set='6-31g*') file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') file.write(' \n') + else: + print(f"Invalid job_name: {job_name}, please use 'cation', 'anion', 'sp_a' or 'sp_c'. See documentation for more details.") \ No newline at end of file diff --git a/build/lib/QCflow/write_xTB.py b/build/lib/QCflow/write_xTB.py new file mode 100644 index 0000000..29c3976 --- /dev/null +++ b/build/lib/QCflow/write_xTB.py @@ -0,0 +1,377 @@ +from QCflow.fragments import * +from QCflow.find_torsion import * +from QCflow.write_gaussian import * +from QCflow.slurm import * +from QCflow.torsion_parser import * +from QCflow.run_gaussian import * +from QCflow.write_psi4 import * +import os +import json +import subprocess +import tempfile + +def run_xtb_optimisation(xyz: str, output_file: str, charge=0, uhf=0, gfn=2): + """ + Runs a xTB geometry optimisation using GFN2-xTB and saves the result to output_file. + + Parameters + ---------- + xyz (str): The input geometry in XYZ format as a string. + output_file (str): The path to the file where the optimised geometry will be saved + charge (int, optional): The total charge of the system. Default is 0. + uhf (int, optional): The number of unpaired electrons (spin multiplicity - 1). Default is 0 (singlet). + gfn (int, optional): The GFN-xTB method to use (1, 2, or 3). Default is 2 (GFN2-xTB). + + Returns + ------- + None: Writes the optimised geometry to output_file. + + """ + with tempfile.TemporaryDirectory() as tmpdir: + xyz_path = os.path.join(tmpdir, "input.xyz") + with open(xyz_path, 'w') as f: + f.write(xyz) + + cmd = [ + "xtb", xyz_path, + "--opt", + "--chrg", str(charge), + "--uhf", str(uhf), + "--gfn", str(gfn) + ] + + try: + subprocess.run(cmd, cwd=tmpdir, check=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) + except subprocess.CalledProcessError as e: + raise RuntimeError(f"xTB failed: {e.stderr.decode()}") + + opt_path = os.path.join(tmpdir, "xtbopt.xyz") + if not os.path.exists(opt_path): + raise FileNotFoundError("Optimised structure (xtbopt.xyz) not found.") + + with open(opt_path, 'r') as f: + optimised_xyz = f.read() + + with open(output_file, 'w') as f: + f.write(optimised_xyz) + + print(f"Optimised structure saved to: {output_file}") + +def load_xyz_from_file(filepath): + """ + Loads an XYZ file and returns the content as a string. + + Parameters + ---------- + filepath (str): The path to the XYZ file to be loaded. + + Returns + ------- + str: The content of the XYZ file as a string. + + """ + with open(filepath, 'r') as f: + lines = f.readlines() + + # Remove the first two lines (atom count and comment/energy) + coordinates = lines[2:] + + # Recalculate number of atoms just in case + num_atoms = len(coordinates) + + # Build the xyz string + xyz_string = f"{num_atoms}\n\n" + "".join(coordinates) + return xyz_string + + +def write_xTB_psi4(job_name, mol_name, smile, functional='b3lyp', basis_set='6-31g*', mol=None, conformer=None): + """ + Runs a xTB geometry optimisation and caulcates electronic proeprities using Psi4, based on the provided parameters. + + Parameters + ---------- + job_name (str): The type of job to run. Possible values: + - 'opt': Optimisation neutral + mol_name (str): The name of the molecule. + smile (str): The SMILE string of the molecule. + functional (str, optional): The functional to use. Default is 'b3lyp'. + basis_set (str, optional): The basis set to use. Default is '6-31g'. + mol (rdkit.Chem.rdchem.Mol, optional): The RDKit embedded molecule object. + conformer (rdkit.Chem.rdchem.Conformer, optional): The RDKit conformer of the molecule. + + Returns + ------- + None: Writes the psi4 input file to disk. + + """ + if (job_name=='opt'): + + old_chk = ' ' + symbols = [a.GetSymbol() for a in mol.GetAtoms()] + # get x y z coords + geometry = conformer + calculation = 'opt' + torsion_data = f' \n' + + if Descriptors.NumRadicalElectrons(mol) == 0: + mult_chg = '0 1' + if Descriptors.NumRadicalElectrons(mol) == 1: + mult_chg = '0 2' + if Descriptors.NumRadicalElectrons(mol) == 2: + mult_chg = '0 3' + + file_name = f'{mol_name}_{job_name}_xTB.py' + #Includes information about basis set to allow for ramping + title = f'{mol_name} {job_name} Smile String: {smile}' + number_proc = 'set_num_threads(4)' + memory_num = '4000 mb' + number_of_atoms = mol.GetNumAtoms() + + with open(file_name, 'w') as file: + file.write('import subprocess\n') + file.write('import tempfile\n') + file.write('import os\n') + file.write('import psi4\n') + file.write('from QCflow.write_xTB import *\n') + file.write(' \n')# + file.write(f'xyz_string = """{number_of_atoms}\n') + file.write(' \n') + if old_chk==' ': + for atom,symbol in enumerate(symbols): + p = geometry.GetAtomPosition(atom) + # atom x y z + line = f' {symbol} {p.x:.5f} {p.y:.5f} {p.z:.5f} \n' + file.write(line) + file.write('""" \n') + file.write(' \n') + file.write('run_xtb_optimisation(xyz_string, "optimised_structure.xyz", gfn=2)\n') + file.write(' \n') + file.write(f"mol_sp = psi4.geometry(open('optimised_structure.xyz').read())\n")#use geometry from opt + file.write('psi4.core.set_active_molecule(mol_sp)\n')#set as active molecule + file.write(f"mol_sp.set_molecular_charge(0) \n")# set charge + file.write(f"mol_sp.set_multiplicity(1) \n")# set multiplicity + file.write(' \n') + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df'}})\n") + file.write(' \n') + file.write(f"energy, wfn = psi4.energy('{functional}', return_wfn=True) \n") + file.write(' \n') + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('hartree_to_ev = 27.2114079527\n') + file.write('energy_ev = energy * hartree_to_ev\n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_xTB_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Single Point energy: {energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + file.write(' \n') + + +def write_xTB_psi4_reorg(job_name, mol_name, functional='b3lyp', basis_set='6-31g*'): + """ + Generates a xTB and psi4 input file based on the provided parameters for reorganisation calculations. + + Parameters + ---------- + job_name (str): The type of job to run. Possible values: + - 'cation': Geometry optimisation cation (opt_c) and single of neutral charge, cation geometry (n_c_geo) + - 'anion': Geometry optimisation anion (opt_a) and single of neutral charge, anion geometry (n_a_geo) + - 'sp_c': Single point calculation of neutral geometry at cation charge + - 'sp_a': Single point calculation of neutral geometry at anion charge + mol_name (str): The name of the molecule. + functional (str, optional): The functional to use. Default is 'b3lyp'. + basis_set (str, optional): The basis set to use. Default is '6-31g'. + + Returns + ------- + None: Writes the psi4 input file to disk. + + """ + if (job_name=='cation'): + + charge = 1 + mult = 2 + + job_type = 'opt_c' + + job_sp = 'n_c_geo' + + with open('optimised_structure.xyz', 'r') as f: + lines = f.readlines() + number_of_atoms = int(lines[0].strip()) + coordinate_lines = [line.strip() for line in lines[2:2 + number_of_atoms]] + + + if (job_name=='anion'): + + charge = -1 + mult = 2 + + job_type = 'opt_a' + + job_sp = 'n_a_geo' + + with open('optimised_structure.xyz', 'r') as f: + lines = f.readlines() + number_of_atoms = int(lines[0].strip()) + coordinate_lines = [line.strip() for line in lines[2:2 + number_of_atoms]] + + if (job_name=='sp_c'): + + charge = 1 + mult = 2 + + with open('optimised_structure.xyz', 'r') as f: + lines = f.readlines() + number_of_atoms = int(lines[0].strip()) + coordinate_lines = [line.strip() for line in lines[2:2 + number_of_atoms]] + + if (job_name=='sp_a'): + + charge = -1 + mult = 2 + + with open('optimised_structure.xyz', 'r') as f: + lines = f.readlines() + number_of_atoms = int(lines[0].strip()) + coordinate_lines = [line.strip() for line in lines[2:2 + number_of_atoms]] + + + file_name = f'{mol_name}_{job_name}.py' + #Includes information about basis set to allow for ramping + memory_num = '4000 mb' + + if (job_name=='anion') or (job_name=='cation'): + + with open(file_name, 'w') as file: + file.write('import subprocess\n') + file.write('import tempfile\n') + file.write('import os\n') + file.write('import psi4\n') + file.write('from QCflow.write_xTB import *\n') + file.write(' \n')# + file.write(f'xyz_string = """{number_of_atoms}\n') + file.write(' \n') + # Now insert the parsed coordinate lines + for line in coordinate_lines: + file.write(f'{line}\n') + file.write('""" \n') + file.write(' \n') + file.write(f"run_xtb_optimisation(xyz_string, '{job_type}_optimised_structure.xyz', charge={charge}, uhf={mult}, gfn=2)\n") + file.write(' \n') + file.write(f"mol = psi4.geometry(open('{job_type}_optimised_structure.xyz').read())\n")#use geometry from opt + file.write('psi4.core.set_active_molecule(mol)\n')#set as active molecule + file.write(f"mol.set_molecular_charge({charge}) \n")# set charge + file.write(f"mol.set_multiplicity({mult}) \n")# set multiplicity + file.write("psi4.set_module_options('scf', {'maxiter': 100}) \n") + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df', 'reference': 'uhf'}})\n") + file.write(' \n') + file.write(f"energy, wfn = psi4.optimize('{functional}', engine='geometric', return_wfn=True) \n") + file.write(' \n') + file.write("optimized_geometry_xyz = wfn.molecule().to_string('xyz') \n") + file.write(f"with open('{mol_name}_{job_type}.xyz', 'w') as xyz_file:\n") + file.write(" xyz_file.write(optimized_geometry_xyz) \n") + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('hartree_to_ev = 27.2114079527 \n') + file.write('energy_ev = energy * hartree_to_ev\n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_{job_type}_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Optimized energy: {energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + file.write(' \n')# Now do the n_c_geo or n_a_geo part + file.write(f"mol_sp = psi4.geometry(open('{mol_name}_{job_type}.xyz').read())\n")#use geometry from opt + file.write('psi4.core.set_active_molecule(mol_sp)\n')#set as active molecule + file.write(f"mol_sp.set_molecular_charge(0) \n")# set charge + file.write(f"mol_sp.set_multiplicity(1) \n")# set multiplicity + file.write(' \n') + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df'}})\n") + file.write(' \n') + file.write(f"energy, wfn = psi4.energy('{functional}', return_wfn=True) \n") + file.write(' \n') + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('energy_ev = energy * hartree_to_ev\n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_{job_sp}_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Single Point energy: {energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + file.write(' \n') + + if (job_name=='sp_c') or (job_name=='sp_a'): + with open(file_name, 'w') as file: + file.write('import psi4\n') + file.write(' \n')# + file.write(f"psi4.set_memory('{memory_num}')\n") + file.write('psi4.set_num_threads(4)\n') + file.write(' \n')# + file.write(f"mol = psi4.geometry(open('optimised_structure.xyz').read())\n")#use geometry from opt + file.write('psi4.core.set_active_molecule(mol)\n')#set as active molecule + file.write(f"mol.set_molecular_charge({charge}) \n")# set charge + file.write(f"mol.set_multiplicity({mult}) \n")# set multiplicity + file.write(' \n') + file.write("psi4.set_module_options('scf', {'maxiter': 50}) \n") + file.write(f"psi4.set_options({{'basis': '{basis_set}', 'scf_type': 'df', 'reference': 'uhf'}})\n") + file.write(' \n') + file.write(f"energy, wfn = psi4.energy('{functional}', return_wfn=True) \n") + file.write(' \n') + file.write(' \n') + file.write(f"psi4.core.set_active_molecule(wfn.molecule()) \n") + file.write(' \n') + file.write('orbital_energies = wfn.epsilon_a_subset("AO", "ALL").np \n') + file.write(' \n') + file.write('n_occ_alpha = wfn.nalpha() \n') + file.write('homo_energy = orbital_energies[n_occ_alpha - 1] \n') + file.write('lumo_energy = orbital_energies[n_occ_alpha] \n') + file.write('energy_gap = lumo_energy - homo_energy \n') + file.write(' \n') + file.write('hartree_to_ev = 27.2114079527 \n') + file.write('energy_ev = energy * hartree_to_ev\n') + file.write('homo_energy_ev = homo_energy * hartree_to_ev \n') + file.write('lumo_energy_ev = lumo_energy * hartree_to_ev \n') + file.write('energy_gap_ev = lumo_energy_ev - homo_energy_ev \n') + file.write(' \n') + file.write(f"with open('{mol_name}_{job_name}_energy_and_gap.txt', 'w') as file:\n") + file.write(' file.write(f"Single Point energy: {energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"HOMO: {homo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"LUMO: {lumo_energy_ev:.6f} eV\\n") \n') + file.write(' file.write(f"Energy gap (HOMO-LUMO): {energy_gap_ev:.6f} eV\\n") \n') + file.write(' \n') diff --git a/docs/api.rst b/docs/api.rst index 05811a5..2c78cec 100755 --- a/docs/api.rst +++ b/docs/api.rst @@ -35,6 +35,15 @@ QCflow.write_psi4 :show-inheritance: :no-imported-members: +QCflow.write_xTB +===================== + +.. automodule:: QCflow.write_xTB + :members: + :undoc-members: + :show-inheritance: + :no-imported-members: + QCflow.slurm ============ @@ -62,6 +71,15 @@ QCflow.run_psi4 :show-inheritance: :no-imported-members: +QCflow.run_xTB +=================== + +.. automodule:: QCflow.run_xTB + :members: + :undoc-members: + :show-inheritance: + :no-imported-members: + QCflow.testing_data =================== diff --git a/docs/conf.py b/docs/conf.py index 5df0397..03cd8da 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -25,7 +25,7 @@ project = 'QCflow' author = 'Tristan Stephens-Jones, Micaela Matta' -release = '0.4.0' +release = '0.5.0' # -- General configuration --------------------------------------------------- diff --git a/qcflow.egg-info/PKG-INFO b/qcflow.egg-info/PKG-INFO index 7eadd1a..1b7fd78 100644 --- a/qcflow.egg-info/PKG-INFO +++ b/qcflow.egg-info/PKG-INFO @@ -1,5 +1,5 @@ Metadata-Version: 2.4 Name: qcflow -Version: 0.4.0 +Version: 0.5.0 License-File: LICENSE Dynamic: license-file diff --git a/qcflow.egg-info/SOURCES.txt b/qcflow.egg-info/SOURCES.txt index b29fdcc..02b3ee5 100755 --- a/qcflow.egg-info/SOURCES.txt +++ b/qcflow.egg-info/SOURCES.txt @@ -8,12 +8,14 @@ QCflow/fragments.py QCflow/load_gaussian.py QCflow/run_gaussian.py QCflow/run_psi4.py +QCflow/run_xTB.py QCflow/sa_score.py QCflow/slurm.py QCflow/testing_data.py QCflow/torsion_parser.py QCflow/write_gaussian.py QCflow/write_psi4.py +QCflow/write_xTB.py qcflow.egg-info/PKG-INFO qcflow.egg-info/SOURCES.txt qcflow.egg-info/dependency_links.txt diff --git a/setup.py b/setup.py index 3d3e24b..277da9b 100755 --- a/setup.py +++ b/setup.py @@ -3,6 +3,6 @@ # this can be more complicated and clever, but for now it's minimal setup( name='qcflow', - version='0.4.0', + version='0.5.0', packages=find_packages(include=['QCflow', 'QCflow.*']), )