New Features of DelPhi Program
Below are described some of the new features that have been introduced to DelPhi:
- Multi-dielectric code
- Geometric objects
- Geometric-shaped charge distributions
- Non-linear PBE solving routine
- Multi-salt code
- Energy partitioning
- Other general issues
The original version of Delphi program was able to handle a scenario where a molecule, e.g. a set of atoms whose radius and charge were picked up from suitable files (.siz and .crg), was immersed in a solution. So, only a "two-media" world was considered.
In the present version, many different objects can be inserted in the scenario. They can be either sets of atoms coming from .pdb file or geometric objects; each of them can have its own dielectric constant.
Special attention has been paid to the compatibility with the old version of DelPhi; therefore one can use the new program in the same way as before. On the other hand, one can also take advantage of the new possibilities and build the desired scenario by writing in the parameter file the statement insobj. Then the user will be guided through the insertion of molecules, objects and charge distributions. A suitable fort.13 will be the result of this process and the following time it will no longer be necessary to do the construction phase, if one wants to keep the scenario as it is.
A complete description of a molecule is not always the most convenient way to represent it; this could be because a detailed calculation is computationally too heavy, or rather because one might prefer a schematic description of a system one doesn't know very well in order to have semi-quantitative results.
These are the main reasons that led to endow the new DelPhi code with the capability of inserting geometrical dielectric objects and charge distributions of simple geometrical shapes.
In this context, the user should be allowed to insert various dielectric objects of different types, being a molecule just one of the possible choices. In order to accomplish this, a set of data structures has been built that contains:
- an ordinal number, that labels the object,
- a field that contains the dielectric constant of the medium constituting the object
- an integer called "objectype" that tells if the object is a molecule (0), a sphere (1), a cylinder (2), a cone (3) or a box (4).
- the parameters related to its dimension and position, which are:
- for a molecule : name of the .pdb file,
- for a sphere: center and radius,
- for a cylinder: centers of the two bases and radius,
- for a cone: center of the basis, vertex, opening angle,
- for a box: four adjacent vertices disposed as the origin and the three vectors of a right-handed vector basis (x, y and z, respectively).
There is also the possibility to create more complex shapes: if two objects turn out to be overlapping, the overlapping region will have the dielectric constant of the object assigned later. Exploiting this fact, one can build unions or intersections of the mentioned objects, each of them having the desired dielectric constant, ending up with a more complex scenario.
As a rule, all molecules superimpose on the other objects regardless of the order in which they have been inserted.
It is under development the possibility of using a different probe radius to build the MS on the part of any molecule that is buried within an object, this value is assigned through the radpr2=[value] statement; by default, it is stated to be equal to the "regular" probe radius.
The insobj keyword, put in the parameter file, tells the program to start a routine that asks for the set of molecules, objects and charge distributions that compose the system.
At present, the boundary conditions are not engineered to comply with a system where any object is lying partly outside the cube where the PBE is being solved. If the user is interested to such a system he can start with a low percentage filling where the usual boundary conditions are sufficient. Then he can slowly increase the perfil until it gets higher than 100% using the focussing method.
Geometric-shaped charge distributions
Charge is a property of matter, but one can in principle consider punctual charges or other distributions, also in the empty space. On the other side, polarization charge induced in the matter by the presence of free charge depends essentially from the dielectric where the free charge is embedded in.
In order to handle both these situations, in the new DelPhi program, various types of uniform charge distributions can be inserted and the user can specify if these belong to an object or not. As already mentioned, this is really useful when we are dealing with a charge very close to an interface between different dielectric media and we want to be sure the program assigns the right dielectric to it.
The charge distributions that have been taken into consideration may have the same overall shapes of the objects plus some other.
In fact, a data structure is built to contain:
- an ordinal number, that labels the distribution,
- a field that tells if the distribution is within the volume or only onto the surface,
- a field that tells if the distribution is "free" or belongs to one of the previously inserted objects,
- a field which contains the total charge amount,
- a number called "distrtype" that tells if the shape is spherical (1), cylindrical (2), conic (3), box-shaped (4), point charge (8), linear (9), disk-shaped (10), rectangular (11),
- the parameters related to its dimension and position, which are the same as corresponding objects if "distrtype" is less or equal to 4, and:
- the point coordinates for number 8,
- the two extremes for number 9,
- the disk center and three points disposed as the three vectors of a right-handed vector basis (x, y and z, respectively), and the disk radius for number 10,
- three adjacent vertices disposed as the origin and the two vectors of a right-handed two dimensional vector basis (x and y, respectively), for number 11.
The same routine activated by the statement insobj takes care of charge distributions.
An algorithm has been implemented to be able to solve the non-linear PBE, as described in J. Phys. Chem. B, 105, 6507-6514(2001).
For a fast and reliable convergence of the non-linear equation a good relaxation parameter has to be used. A heuristic algorithm is used to estimate it; nonetheless, for some particularly pathological condition, it is possible that the choice is not optimal and that the convergence is very slow or even absent.
In order to deal also with those systems, the new statement relpar=[value] can be put in the parameter file to manually assign the relaxation parameter.
A general set of guidelines to the choice of the relaxation parameter follows:
- The optimal relaxation parameter for non-linear system generally, not always, turns out to be less than 1. Pathological conditions, like a very high charge in proximity of the surface, ionic strength close to zero but not null and some particular scale values, can require a parameter much lower, down to 0.001 or even lower. On the other hand, a too low value has the drawback to slow down the convergence and, in the worst case, to give only an appearance of convergence. Thus, the advice is to start with a value not too low and then to decrease it until the convergence is reached.
The valence of the salt that is in solution can really make a difference for the electric field in solution, influencing the field also at some distance from the molecule.
Taking into account this phenomenon requires handling the more general expression for the charge distribution in solution, to consider all the possible terms in it, and finally to substitute them in the original PBE.
Some new statements have been added in order to handle this new situation:
- salt=[conc first kind of salt, moles/liter], val+1=[valence of positive ion in salt type 1], val-1=[valence of negative ion in salt type 1];
- salt2=[conc second kind of salt, moles/liter], val+2=[valence of positive ion in salt type 2], val-2=[valence of negative ion in salt type 2];
A new partition of the energy is considered, as described J. Phys. Chem. B, 105, 6507-6514(2001). Accordingly, the electrostatic energy is subdivided in: coulombic, reaction field, self-reaction field, external ion contribution, osmotic term, "rho*phi". The three last term appear only when the ionic strength is different from zero and the last two terms cancel out in the linear approximation.
Coulombic and reaction field energies are calculated as in the previous code.
The osmotic and "rho*phi" terms are calculated whenever the program is asked to run nonlinear iterations.
The calculation of the external ion contribution is triggered by the flag ion in the energy statement; it calculates the direct interaction of ions to real charges.
This direct calculation has some intrinsic drawbacks: first of all, it is computationally very expensive, secondly, this energy is lacking of all the ionic contribution relative to ions located outside the box. Another way to get an estimation of the same energy term is to subtract the grid and reaction field energy in two cases: with and without salt. In case of non-linear PBE also the osmotic and electrostatic stress ("rho*phi") have to be subtracted.
The self-reaction field energy is calculated whenever the usual reaction field energy is calculated; it is an attempt to estimate the energy that the system gains when a point charge is moved from vacuum to a dielectric region. In this model, it is assumed that polarization charge builds up with spherical symmetry at a distance the user can decide via the statement radpolext=[radius], which by default is set to be 1Å. This information is no attempt to give a quantitative estimation of this energy but just to take into account in a more detailed way the fact that the most favorable place for a charge is within a high dielectric medium. This necessity arises when a system with many different dielectric regions is studied and it is no longer straightforward to decide what is the reference state.
Other general issues
- Two more statements have been added: rmsc=[conc. threshold moles/liter] and maxc= conc. threshold moles/liter] that put two new convergence criteria. The first concerns the rms change and the second the maximum change. They can be assigned also together with the number of iterations; the program will stop iterating as soon as one of the criteria will be met.
- Containment of memory requirement for large system, being able to run membrane with about 50000 atoms to grid size 395 on Octane 12000 with 2 Gb Ram (32 bit) and up roughly to 571 grid size for a 64 bit version.
- Porting of the code to Linux operating system, compiled using Portland Group compiler.
- Porting of the code onto Windows Operating System, using Microsoft Developer Studio compiler.
- Porting of the code on SGI, IRIX 64 bit.