Skip to content

Commit

Permalink
Fixed RESP for 2-atom molecules
Browse files Browse the repository at this point in the history
  • Loading branch information
tkpiskorz committed Aug 30, 2024
1 parent 1d61f4c commit 7b60482
Show file tree
Hide file tree
Showing 8 changed files with 72 additions and 42 deletions.
20 changes: 15 additions & 5 deletions docs/source/examples.rst
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,11 @@ Examples
This is extended desription of the examples provided in metallicious package (see `examples <https://github.com/tkpiskorz/metallicious/tree/main/metallicious/examples>`_)


Example 1: Quick start (see `example_1/ <https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/example1_quick_start>`_)
Example 1: Quick start
------------

(see `example_1/ <https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/example1_quick_start>`_)

Simple parametrization of the cage with two metals. As input, we provide a coordination file in PDB format and non-bonded
force-field parameters in *.top format (GROMACS). Moreover, LJ parameters of the Ru and Pd are taken from the universal force-field
(vdw_type='uff'), as Ru is available only in this library.
Expand All @@ -21,9 +23,12 @@ force-field parameters in *.top format (GROMACS). Moreover, LJ parameters of the
cage.parametrize(out_coord='out.pdb', out_topol='out.top')
Example 2a: Automated initial force-field parameters (see `example_2/ <https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/example2_no_topology>`_) If the force-field parameters are missing, the linker(s) can be parametrized using AmberTools with a GAFF force field.
Example 2a: Automated initial force-field parameters
------------

(see `example_2/ <https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/example2_no_topology>`_)

If the force-field parameters are missing, the linker(s) can be parametrized using AmberTools with a GAFF force field.
*metallicious* will prepare the non-bonded force-field parameters of the whole structure, which are saved to the folder "init_topol."

.. code-block:: python
Expand Down Expand Up @@ -55,9 +60,11 @@ Alternatively, you prepare_initial_topology can be specified in .parametrize:
cage.parametrize(out_coord='out.inpcrd', out_topol='out.prmtop', prepare_initial_topology=True)
Example 3: Homoleptic cage (see `example_3/ <https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/example3_only_linker_topology>`_)
Example 3: Homoleptic cage
------------

(see `example_3/ <https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/example3_only_linker_topology>`_)

In the case of homoleptic cages (all linkers are the same), *metallicious* can use force-field parameters of single linker,
which will be used to generate the initial topology of the whole structure:

Expand Down Expand Up @@ -97,9 +104,11 @@ Two solutions are available:
1. Create a new template, which is automated but time-consuming
2. use truncation schemes, which are fast but with reduced accuracy (caution is also needed)

Example 4: Parametrization of new template (see `example_4/ <https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/example4_template_parametrization>`_)
Example 4: Parametrization of new template
------------

(see `example_4/ <https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/example4_template_parametrization>`_)

If template parameters are unavailable, one may want to parametrize them. In *metallicious*, this is done by
specifying the multiplicity of the metal, which signals *metallicious* to perform QM calculations if necessary. For this functionality, the additional
dependencies (see installation guide) are needed (`autode <https://github.com/duartegroup/autodE>`_, `ORCA <https://orcaforum.kofo.mpg.de/app.php/portal>`_, and `psiRESP <https://github.com/lilyminium/psiresp>`_).
Expand Down Expand Up @@ -134,9 +143,10 @@ which can be changed by specifying "keywords" in the supramolecular_structure cl
cage.parametrize(out_coord='out.pdb', out_topol='out.top', prepare_initial_topology=True)
Example 5: Truncation schemes (see `example_5/ https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/example5_truncation_scheme`_)
Example 5: Truncation schemes
------------

(see `example_5/ https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/example5_truncation_scheme`_)

Truncation schemes allow the "recycling" of existing templates from the library by reducing their size, which might match the metal site of interest.
Three schemes are available, cutting the template at a distance of 3-bond, 2-bond, and 1-bond from the metal centre.
Expand Down
87 changes: 53 additions & 34 deletions docs/source/tutorials.rst
Original file line number Diff line number Diff line change
Expand Up @@ -95,76 +95,95 @@ This requires only change of the library_directory, as there are no templates in
Tutorial 5: GROMACS tutorial for MD simulations of metallorganic structures from scratch
-----------


This tutorial will guide you how to perform whole MD simulations from scratch. We will simulate a Ga4L6 cage in water.
Firstly, let's do the topology of the cage from *.xyz structure:
This tutorial will guide you in performing an MD simulation from scratch. We will simulate a Fujita's Pd4L6 cage in water.
The files needed for this tutorial are available `here <https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/tutorial5/tutorial.zip>`_.
Firstly, let's do the topology of the cage from *.xyz structure using metallicious:
.. code-block:: python
from metallicious import supramolecular_structure
cage = supramolecular_structure('cage_start.pdb', metal_charges={'Pd': 2}, LJ_type='uff')
cage.parametrize(out_coord='0_cage.pdb', out_topol='topol.top', prepare_initial_topology=True)
Execule the program:
Execute the program:

.. code-block:: bash
python script.py
With the topology read, let's setup the system. Firstly let's modify the size of the simulation box:
With the topology read, we will set up the system. Firstly, let's modify the size of the simulation box to accommodate
the cage (we use 2 nm from the cage to the wall of the simulation box):

.. code-block:: bash
gmx editconf -f 0_cage.pdb -d 2 -o 1_box.pdb
We solvate the system with water (we usee 3-site water model):
gmx editconf -f 0_cage.pdb -d 2 -o 1_box.pdb
We solvate the system with water (we will use a 3-site water model):

.. code-block:: bash
gmx solvate -cp 1_box.pdb -cs spc216.gro -o 2_solv.gro -p topol.top
gmx solvate -cp 1_box.pdb -cs spc216.gro -o 2_solv.gro -p topol.top
Now, the topology file needs to be changed. Open the "topol.top" file. Firstly, we will change the header:

.. code-block:: bash
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 yes 0.5 0.833333
to:
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 yes 0.5 0.833333
#include "amber99sb-ildn.ff/forcefield.itp"
to force-field parameters, which include atom types of water:

.. code-block:: bash
#include "amber99sb-ildn.ff/forcefield.itp"
This will include all LJ parameters for atoms (including solvent and ions). Next, add information about the force-field
parameters for water and ions just under the section "[ atomtypes ]" (just
before the first [ moleculetype ]). In this case, we will use SPC/E water model:

.. code-block:: bash
This will include all LJ paramters for atoms (including solvent and ions). Under the section "[ atomtypes ]" (just
before first [ moleculetype ]) we add information about the force-field paramters for water and ions topology.
In this case we will use SPC/E water model:
#include "amber99sb-ildn.ff/spce.itp"
#include "amber99sb-ildn.ff/ions.itp"
#include "amber99sb-ildn.ff/spce.itp"
#include "amber99sb-ildn.ff/ions.itp"
Now, we need to add ions. Firstly, we need to create a *.tpr file:
Now we need to add ions. Firstly, we need to create a tpr file
.. code-block:: bash
gmx grompp -f em.mdp -c 2_solv.gro -p topol.top -o ions.tpr -maxwarn 2
gmx grompp -f em.mdp -c 2_solv.gro -p topol.top -o ions.tpr -maxwarn 2
Then, let's add ions:

Then you add ions:
.. code-block:: bash
gmx genion -s ions.tpr -o 3_ions.gro -p topol.top -pname NA -nname CL -neutral
gmx genion -s ions.tpr -o 3_ions.gro -p topol.top -pname NA -nname CL -neutral
Select "SOL" to replace bulk water with the counterions.
We minimize the created system (to "remove bad contacts") by firstly creating *tpr file, and then performing miniization.
We minimize the created system (to "remove bad contacts") by making a *.tpr file and then performing minimization.
gmx grompp -f em.mdp -c 3_ions.gro -p topol.top -o 4_em.tpr
gmx mdrun -v -deffnm 4_em
.. code-block:: bash
Equilibration using NVT and NPT:
gmx grompp -f em.mdp -c 3_ions.gro -p topol.top -o 4_em.tpr
gmx mdrun -v -deffnm 4_em
gmx grompp -f nvt.mdp -c 4_em.gro -p topol.top -o 5_nvt.tpr
gmx mdrun -v -deffnm 5_nvt
We equilibrate system using NVT and NPT ensemble:

.. code-block:: bash
gmx grompp -f npt.mdp -c 5_nvt.gro -p topol.top -o 6_npt.tpr
gmx mdrun -v -deffnm 6_npt
gmx grompp -f nvt.mdp -c 4_em.gro -p topol.top -o 5_nvt.tpr
gmx mdrun -v -deffnm 5_nvt
gmx grompp -f npt.mdp -c 5_nvt.gro -p topol.top -o 6_npt.tpr
gmx mdrun -v -deffnm 6_npt
Finally, we perform production run (ideally done on a performance cluster (HPC) with GPUs):

.. code-block:: bash
And finally production run (ideally done on high performance cluster (HPC) with GPUs):
gmx grompp -f run.mdp -c 6_npt.gro -p topol.top -o 7_run.tpr
gmx mdrun -v -deffnm 7_run
gmx grompp -f run.mdp -c 6_npt.gro -p topol.top -o 7_run.tpr
gmx mdrun -v -deffnm 7_run
As result you obtaine final frame 7_run.gro and trajectory 7_run.xtc, which you can visualise using for example VMD.
As a result, we should obtain the final frame 7_run.gro and trajectory 7_run.xtc, which you can visualise using a
molecular visualization program (i.g., VMD). For your convenience, `here <https://github.com/duartegroup/metallicious/tree/main/metallicious/examples/tutorial5/final.zip>`_ are the final files compliled.
2 changes: 1 addition & 1 deletion metallicious.egg-info/PKG-INFO
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Metadata-Version: 2.1
Name: metallicious
Version: 0.2.32
Version: 0.2.33
License-File: LICENSE.md
Requires-Dist: numpy
Requires-Dist: MDAnalysis
Expand Down
2 changes: 2 additions & 0 deletions metallicious.egg-info/SOURCES.txt
Original file line number Diff line number Diff line change
Expand Up @@ -114,6 +114,8 @@ metallicious/library/Fe_2_merz-opc_3.pdb
metallicious/library/Fe_2_merz-opc_3.top
metallicious/library/Fe_2_uff_0.pdb
metallicious/library/Fe_2_uff_0.top
metallicious/library/Fe_2_uff_1.pdb
metallicious/library/Fe_2_uff_1.top
metallicious/library/Fe_2_zhang-opc_0.pdb
metallicious/library/Fe_2_zhang-opc_0.top
metallicious/library/Fe_2_zhang-opc_1.pdb
Expand Down
Binary file added metallicious/examples/tutorial5/final.zip
Binary file not shown.
Binary file added metallicious/examples/tutorial5/tutorial.zip
Binary file not shown.
1 change: 0 additions & 1 deletion metallicious/parametrize_new_sites.py
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,6 @@ def __init__(self, filename, metal_charge_mult=None, metal_charges=None, LJ_type

elif LJ_type is None:
for LJ_type in vdw_data:

present = [(metal_name in vdw_data[LJ_type] or f"{metal_name:}{self.metal_charge_dict[metal_name]:}" in vdw_data[LJ_type]) for metal_name in self.metal_names]
if sum(present) == len(present):
self.vdw_type = LJ_type
Expand Down
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"

[project]
name = "metallicious"
version = "0.2.32"
version = "0.2.33"

dependencies = [
"numpy",
Expand Down

0 comments on commit 7b60482

Please sign in to comment.