Generating Free Energy Profiles ================================= This tutorial will guide you through the process of generating energy profiles as a collective variable (CV) is modulated with a calculation. Learning Objectives ------------------- - Learn how to generate energy profiles with different electronic structure interfaces. - Learn how to generate energy profiles from DataBase queries. Required Files -------------- - The database must be previously generated, and queryable. - The tutorial python script is located in examples/EnergyProfiles Tutorial -------- To generate the energy profiles, we will use the `calculate_profiles.py` script. This is a bit more complex than the previous example, but it is still fairly short. To start, we will import the necessary modules. .. literalinclude:: ../../../../examples/EnergyProfiles/calculate_profiles.py :language: python :end-before: justplot = Now, this script has two modes to be run because it can be a bit slow to run. To start, run the script with .. code-block:: python justplot = False With this option, it will actually run calculations. Note - if it can't find the xtb or psi4 scans, it will set this back to false and run the calculations so, by default, you can run it with justplot = True and it will run the calculations if they are not found. Ethane-Like Torsional Scan ^^^^^^^^^^^^^^^^^^^^^^^^^^ We will start with a small, ethane-like molecule with a torsional scan. .. literalinclude:: ../../../../examples/EnergyProfiles/calculate_profiles.py :language: python :start-after: if not justplot :end-before: Now lets try it again In this code, it generates an ASE Atoms object (with fake coordinates). Then the newscan object calls the XTB Interface and scans the dihedral angle between atoms 1 2 3 4 from 0 degrees to 360 degrees in 20 degree increments. Finally, ethane_scan.xyz is written and can be visualized with VMD. Torsional Scanning with XTB and a Database query ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Now, we will try generating a scan from the database query instead. We will start by loading the database, generating a query, and then drawing the molecule with rdkit. .. literalinclude:: ../../../../examples/EnergyProfiles/calculate_profiles.py :language: python :start-at: Now lets try it again :end-at: img.show() .. image:: figures/torsionalscan_molecule.png :width: 600px :align: center :alt: The molecule used for the torsional scan, this is 2-[(thiophen-2-ylmethyltrisulfanyl)methyl]thiophene. The molecule is a bit more complex than our first example, but it is still a small(ish) molecule. Its IUPAC name is 2-[(thiophen-2-ylmethyltrisulfanyl)methyl]thiophene, and we will be using it to generate a torsional scan where we rotate one of the dihedral torsions in the chain linking the two thiophene rings. We will first start by taking one of the first matches from the query and converting it to an ASE atoms object. .. literalinclude:: ../../../../examples/EnergyProfiles/calculate_profiles.py :language: python :start-at: Convert the molecule :end-at: print(f"Created Atoms object: {atoms}") Now, like before, we will create a newscan object and run the scan, first with XTB. The XTB progrogram is relatively simple to run with the interface built into `pharmaforge`. In this example, we use the GFN-xTB method. .. literalinclude:: ../../../../examples/EnergyProfiles/calculate_profiles.py :language: python :start-after: XTB torsional scan :end-before: End XTB scan When you run the scan, it will show a plot on the screen (if you are normally able to view matplotlib plots). This is the energy profile of the scan. Repeating the scan with Psi4 (Scan Parallelization) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Now, we will do the same thing but with Psi4. The Psi4 interface is similar, but will actually run a DFT calculation. Here we use B3LYP/6-31G, but you can change this to any functional and basis set supported by Psi4. Psi4 supports two methods of parallelization, you can define the number of threads for a single calculation to use with the `num_threads` argument to `Psi4Interface`, or you can use the `num_workers` argument to run multiple calculations in parallel. Timings can be found in the table below. .. literalinclude:: ../../../../examples/EnergyProfiles/calculate_profiles.py :language: python :start-after: Psi4 torsional scan :end-before: End Psi4 scan .. warning:: Some interfaces (for instance XTB) utilizes all the vailable threads on the system, so you should not use the `num_workers` argument with them, by default this will be set to 1. The code has built in checks, and will revert to the serial interface for interfaces that do not support parallelization. .. list-table:: Timings for Parallel Scans :widths: 50 50 50 50 50 50 :align: center :header-rows: 1 * - Program - 1 thread, 1 worker (s) - 4 threads, 1 workers (s) - 4 threads, 3 workers (s) - 1 thread, 10 workers (s) - 1 thread, 12 workers (s) * - XTB - 18.01 - N/A - N/A - N/A - N/A * - Psi4 - 1225.97 - 413.84 - 175.09 - 160.91 - 169.02 Lastly, we save the data so that we can plot it later without rerunning the calculations. .. literalinclude:: ../../../../examples/EnergyProfiles/calculate_profiles.py :language: python :start-at: xtb_data = np.savetxt :end-at: psi_data = np.savetxt Plotting the Final Scan Data ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Now, to plot the data set, we can run the same script with the justplot variable set to True. .. code-block:: python justplot = True And the following code will plot the data. .. literalinclude:: ../../../../examples/EnergyProfiles/calculate_profiles.py :language: python :start-at: np.loadtxt This plot, which you should have produced with the above command, looks like this: .. image:: figures/torsionalscan.png :width: 1200px :align: center :alt: Comparison of XTB and Psi4 energy profiles Full Code --------- .. literalinclude :: ../../../../examples/EnergyProfiles/calculate_profiles.py :language: python