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:
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 (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.
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.
# 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.