Considering the effects of a solvent in interactions

[hadi9827]

Hello,
I am using Quantum-ATK to investigate the adsorption of a protein molecule (Albumin) on a surface of a slab like ZnO or diamond. Because of the large size of the system, I am using forcefield simulations. I did some simulations with the assumption of the protein molecule being in vaccum and over the surface of the slab. Now, I need to consider the protein being in a solvent (not vacuum)  and over the surface of the slab. I intend to go for implicit solvent strategy (although it is not accurate). These are my questions:

1. how I can introduce a solvent to be around my protein molecule instead of the vacuum?

2. Is using the different dielectric constant for the space (currently vacuum) above the surface a good idea for doing that? If yes how can I do that?

3. I used forcefield simulations so far (for the molecule in vacuum). What kind of calculations should I use for the molecule in solvent? Is using forcefield still possible? Generally would you please guide me on choosing the proper calculator in Script Generator for the new case? I need the
calculations be as fast as forcefield if it is possible.

Thank you very much for your responses in advance.

[BradWells]

Hi Hadi,

Thanks for your questions.  In response to your queries:

1)  The simplest way to add a solvent in your system is by using the 'Packmol' tool in the builder.  To use this tool you simply need to define the bulk surface/protein configuration that is to be filled with solvent, and the solvent itself as a molecule configuration.  In the Packmol tool you can then add the bulk and the molecules to be packed in, and the tool will pack the solvent molecules around the surface/protein.  You then may need to then run a short molecular dynamics calculation to fix the periodicity of the solvent, and then you will be ready to start your calculation.

2)  Support for implicit solvent models is a little limited in QuantumATK.  It is possible to use a screened coulomb potential which can mimic some of the effects of implicit solvent (see https://docs.quantumatk.com/manual/Types/CoulombDebye/CoulombDebye.html).  We don't however have a simple solution for implementing models such as Generalized Born Implicit Solvent. 

3)  It is a little hard to give a precise answer here, because you haven't mentioned what property of the system you are trying to calculate.  In general I would recommend using a Forcefield calculator with explicit solvent molecules.  Likely one of the important details of the system will be the differences in the affinity of the protein with the surface compared to solvent.  This can include both entropic and enthalpic components which can only really be modeled with explicit solvent.  The forcefield calculator is also generally the fastest of the available calculators.

If you are able to give some more details on the properties that you are hoping to estimate, I can try to give you a bit more specific advice on how to calculate what you are interested in.

Brad.

[hadi9827]

Hi Brad,
Thank you very much for your very helpful and comprehensive response. Please find the attachment. Thanks.
Hadi

[BradWells]

Hi Hadi, To answer your questions: 2)  I have been looking into setting the dielectric constant in forcefield simulations, and there is a way to directly set the dielectric constant for Coulomb solvers.  In your script, if you have a potential set defined with the variable "potential_set", you can change the dielectric constant with the lines:

Code

dielectric = 1.0
potential_set.addOption(CoulombOption(dielectric / (4*numpy.pi*vacuum_permitivity)))

Here you can change the dielectric to the value that you want.  Unfortunately however, this does not apply to the ReaxFF potentials.  The ReaxFF potential does not use a traditional Coulomb 1/r term, but rather a screened Coulomb expression.  As such, it uses special Coulomb methods that are not altered by this option.  Using ReaxFF, there is no way to change the dielectric constant in the simulation.  In general, adsorbing molecules onto metal or metal-oxide surfaces is a difficult problem, as there are not many potentials that are well parameterized for general metal-organic interactions.  For example, if you look up the original reference paper for the ReaxFF potential that you are using, you will see that the potential has been mostly fitted for Zn-O-Ag interactions.  These parameters have been simply added to existing organic parameters without really considering their interaction.  To get reliable results I think you need to be careful in selecting a potential that is parameterized for the system that you are interested in studying.  In general, it would be much easier simulate adsorption onto a diamond (carbon) surface.  That would eliminate a lot of the problems with finding metal potentials. Another thing to consider is how important modeling actual chemisorption is to your problem.  ReaxFF and other reactive potentials are able to simulate bond breaking and forming, at the cost of a significant slowdown in the simulation.  If you are not as interested in chemisorption, you may be able to simulate the problem using a traditional bonded potential, such as DREIDING.  Bonded potentials are much faster, but cannot simulate bond formation, and therefore cannot simulate chemisorption.  In QuantumATK there is support for DREIDING, but only through the PolymerBuilder and scripting.  DREIDING also contains somewhat appropriate potentials for ZnO and is able to model adsorption of molecules using simple Lennard-Jones and Coulomb functions.  If you think this may suit your calculation I can show you how to set that up in a script.  It is a little tricky combining it with setting the dielectric constant, but it can be done. You can also see documentation for the DREIDING potential at: https://docs.quantumatk.com/manual/Types/DreidingPotentialBuilder/DreidingPotentialBuilder.html 4)  In the tutorials that you mention, the dielectric constant is being set in both DFTB and DFT-LCAO methods.  These calculators perform semi-empirical and full density functional theory calculations respectively.  As such, they are much more computationally intensive than a simple forcefield, as these methods take into account the electronic structure of the material.  Given the size of the albumin molecule and the surface that you need to adsorb it, I would estimate that LCAO calculations are probably too computationally intensive.  It may be possible to use a DFTB calculator with your system, but only to do single point energy calculations and perhaps optimizations of an adsorption position.  It would probably be too slow to run any meaningful MD simulation. Hope that helps you, Brad.

[hadi9827]

Thank you very much for your complete and helpful response.

Unfortunately ReaxFF is the only predefined potential that is available in ATK for my system. I was thinking of defining a potential set for my system manually but I suppose it will be so difficult to do that accurately for getting reliable results.

I will be so grateful if you help me set DREIDING step by step.

In terms of DFTB calculations, I think as you mentioned correctly it is too computationally costly.

Best regards,
Hadi

[BradWells]

Hi Hadi, You can set up the DREIDING potential for your system using the code below.  The easiest thing to do is set up your calculation with the ReaxFF potential, and then in script editor replace the ReaxFF part with the following code.

Code

# Find the bonds in the bulk configuration
bulk_configuration.findBonds()

# Create the potential builder
potential_builder = DreidingPotentialBuilder()

# Assign atom types onto the configuration
potential_builder.assignAtomTypes(bulk_configuration)

# Create a calculator
calculator = potential_builder.createCalculator(bulk_configuration)

# Add a dielectric constant
dielectric = 10
calculator._potential_parameters.addOption(CoulombOption(1/(4*numpy.pi*dielectric*vacuum_permitivity)))

# Set the calculator on the configuration
bulk_configuration.setCalculator(calculator)

# Save the configuration to a file in case it is needed later
nlsave('Input_Configuration.hdf5', bulk_configuration)

To use the DREIDING potential you have to put bonds and atom type tags onto your configuration.  This is handled automatically in the code.  You can also change the value of the 'dielectic' variable to apply the dielectric constant you want. One thing I should mention is that DREIDING uses atomic charges calculated with QEq.  Calculating the charges can take a significant amount of time, depending on the size of your system.  Once you calculate the charges, they are added to the calculator which is saved with the configuration.  If you have to repeat the calculation, you should re-use the same calculator so you don't have to calculate the charges again. Hope this helps, Brad.

[hadi9827]

Thank you very much Brad, for your great responses.