next up previous contents index
Next: To access hidden features Up: eFF Programmers' Guide Previous: To modify current features   Contents   Index


To add new features

Add a new parameter

To add a new parameter, first make an entry for it in the params data structure in eff_global:

  typedef struct
    double new_param;
Then create a parsing function for it in eff_params:
  void ParamNewParam(int argc, char **argv)
    if (argc != 3) error("Usage: new_param = [new parameter]\n");
In the above example, the parameter appears in the input .cfg file as ``new_param = [value]''. However, variations are possible - there may be more arguments, the arguments may be mixtures of strings and integers - which require a more complex parsing function.

Next, associate the parsing function with a parameter name in AddParameters() in eff_params:

  void InitializeParameters()
    AddParameter("new_param", ParamNewParam);

Finally, place an entry in OutputFixedHeader() in eff_output so that the parameter is printed out at the start of the .eff output file and .cfg.restart restart file:

void OutputParams(char *tag, FILE *out_fp)
  fprintf(out_fp, "%s\tnew_param = %f\n", tag, params.new_param);

Now the parameter can be accessed as params.new_param in any file that includes eff_global.h.

Add a new restraint

To add a new restraint, first make a new data structure for it in eff_restraints:
  typedef struct
    int idx1, idx2, idx3;
    double value;
  } RES_NEW;
Here we have assumed that the restraint operates between three particles. If it operates between a variable number of particles, use a structure of the following form instead:
  typedef struct
    int *indices;
    int numindices;
    double value;
  } RES_NEW;

Next, declare a new global variable of type LIST_PTR * to serve as a restraint list, and add statements to initialize and clear it in InitializeRestraints() and ClearRestraints():

  LIST_PTR *restraint_new;             // in the global scope
  restraint_new = new_list_ptr();      // in InitializeRestraints()
  clear_list_ptr(restraint_new);       // in ClearRestraints()

Then create a function for adding the restraint to the restraint list:

  void AddNewRestraint(int idx1, int idx2, int idx3, double value)
    RES_NEW *ptr = (RES_NEW *) malloc(sizeof(RES_NEW));
    ptr->idx1 = idx1;
    ptr->idx2 = idx2;
    ptr->idx3 = idx3;
    ptr->value = value;
    add_list_ptr(restraint_new, ptr);
If the restraint operates on a variable of particles, use AddLineRestraint() as a model for your function instead.

Now we are ready to set up the function which actually enforces the restraint through use of a restraining potential, usually harmonic. Create a function called ApplyNewRestraint():

  double ApplyNewRestraint(int idx1, int idx2, int idx3,
                           double value0, double *value)
     // Gets the positions of the particles
     double x1, y1, z1, x2, y2, z2, x3, y3, z3;
     GetNucElecPosition(idx1, &x1, &y1, &z1);
     GetNucElecPosition(idx2, &x2, &y2, &z2);
     GetNucElecPosition(idx3, &x3, &y3, &z3);

     // Computes the forces on the particles
     double fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3;

     if (value) *value = restraint_measure;

     // Applies forces to the particles
     AddNucElecForce(idx1, fx1, fy1, fz1);
     AddNucElecForce(idx2, fx2, fy2, fz2);
     AddNucElecForce(idx3, fx3, fy3, fz3);

Next add the above function to the main ApplyRestraints() function, making sure it is inserted above the ApplyTetherRestraints() function, which needs to be applied last to work.

  for (ptr = restraint_new->head; ptr != 0; ptr = ptr->next)
    RES_NEW *element = (RES_NEW *) ptr->element;
    restraints_energy += ApplyNewRestraint(
      element->idx1, element->idx2, element->idx3,
      element->value, 0);

Finally, add a similar loop to the OutputRestartRestraints() function:

  for (ptr = restraint_new->head; ptr != 0; ptr = ptr->next)
    RES_NEW *element = (RES_NEW *) ptr->element;
    NucElecString(str1, element->idx1);
    NucElecString(str2, element->idx2);
    NucElecString(str3, element->idx3);
    fprintf(restart_fp, "new %s %s %s %f\n", str1, str2, str3, element->value);

Add a new property to calculate

To calculate a new property, add a function to eff_properties. Then use functions defined in eff_access to get nuclear and electron positions and forces, or functions in eff_dynamics to get velocities and masses.

If the property is only to be computed once, at the end of the calculation, add the property function to OutputEnergy() in eff_output.

Otherwise, if the property is to be computed and printed out every print_every iterations, create a new function in eff_output that outputs the property:

  void OutputNewProperty(FILE *out_fp)
    fprintf(out_fp, "[new_property]\t%f\n", CalcNewProperty();

Use the instructions in Section 4.1 to add a new parameter called ``output_new_property'', added in eff_global as type ``enum OutputType''. Then the user will be able to set the parameter to NONE, ALL, or END to determine how often the property will be printed.

Do a search for OutputPositions in eff_driver, and insert lines of the form:

  if (params.output_new_property == ALL) OutputNewProperty(out_fp);
into the file following the example of OutputPositions, which appears in RunSinglePoint(), MinimizeCallback(), RunMinimize(), and RunDynamics().

Add more degrees of freedom

To add degrees of freedom to the electrons, start by putting extra parameters into the NUCLEUS and ELECTRON typedefs in eff_global.

Next, change the routines ElecNucElec(), ElecElecElec(), and PauliElecElec() in eff_update. There will need to be a way to keep the wavefunction normalized as well.

If the intent is to represent the electrons as linear combinations of primitive functions, it is straightforward enough to expand the kinetic energy and electrostatic sums into terms involving the primitives. An Overlap() function could be defined to aid normalization.

For the Pauli function, one needs to settle on a balance between implicit and explicit orthogonalization. The current functional works via implicit orthogonalization - it makes no attempt to generate orthogonal wavefunctions, but instead estimates the kinetic energy increase that would be expected upon orthogonalization.

With more degrees of freedom, however, there is possibility of changing the functions explicitly so that the wavefunctions are orthogonal. We have had some success using the SNOPT constrained optimizer by Gill, Murray, and Saunders for this purpose. One must be careful not to double-count the orthogonalization energy in that case. The exchange energy term also then becomes more prominent, and may be useful to calculate.

Too much variational flexibility can also lead to linear dependence issues. Two $ s$ floating gaussians can generate a $ p$ like function in the limit that they are near each other - however, this can lead to other numerical problems. This problem can be avoided somewhat by setting appropriate constraints, for example, mandating that primitives have exponents that are in fixed ratios to each other.

The semiclassical equations of motion in eff_dynamics will need to be updated as well, for example with the time dependent variational approach proposed by E. J. Heller, J. Chem. Phys. 64(1), 63-73 (1976).

In eff_minimize, the functions GetMinimizeForces(), GetInitialPosition(), and SetMinimizePositions() provide a mapping from the minimization arrays to the eFF system, and should be updated. In eff_dynamics, the functions GetDynamicsForces(), GetDynamicsPositions(), SetDyanamicsPositions() serve a similar role and should be updated as well.

Our experience

We experimented with using sums of Gaussians to make eFF more accurate. We found that the potential energy surface of $ \mathrm{H_{2}}$ could be made to approach the Hartree-Fock potential surface within a few kcal/mol, using only two floating Gaussian primitives per electron.

However, for systems with more than two electrons, double-counting in the Pauli potential and linear dependencies became problematic. Although there has been success in the literature with using multiple floating Gaussians to represent lone pairs, e.g. A. H. Pakari and F. M. Khalesifard, J. Mol. Struct. Theochem 340:175-183 (1995), those approaches have relied on ``hand-positioned'' Gaussian functions with very limited degrees of freedom. Once these systems were allowed to relax, we observed unphysical geometries and bond energies, as well as linear-dependence issues.

In general, we found it more effective to keep the simple single Gaussian description of electrons, and to add physically motivated terms to the eFF expression to account for other effects, such as the $ p$ character of electrons at a certain range from the nucleus, or electron correlation between opposite spin electrons in close proximity to each other. We have had a good deal of success with this approach, and we believe we have only begun to explore the potential of where these developments could lead.

By keeping the electrons simple and accounting for more complex effects in an implicit way, the equations of motion remain simple, the storage requirements remain light, and eFF remains fast.

Make nuclei floating Gaussians

In the current eFF scheme, nuclei are represented as classical particles. However, protons are light enough that they have a significant quantum delocalization, which can affect proton transfer rates, as in the kinetic isotope effect, and also affect excited electron motions coupled to proton motions.

In the nuclear-electron orbital method proposed by S. P. Webb, T. Iordanov, and S. Hammes-Schiffer, J. Chem. Phys. 117(9), 4106-4118 (2002), hydrogen nuclei as well as electrons are represented by sums of Gaussian functions. The analogous approach in eFF would be to make nuclei as well as electrons floating Gaussian functions.

One way to do this would be to extend the definition of electrons to include non-classical nuclei. These particles would have a kinetic energy inversely proportional to their mass, and would not interact via the Pauli potential (they could, but it would only come about at very small nucleus-nucleus separations).

First, enter a static mass m into the ELECTRON typedef in eff_global. Then modify the function UpdateKineticEnergy() in eff_update, changing the kinetic energy from

$\displaystyle E_{\mathrm{ke}} = \frac{3}{2} \cdot \frac{1}{r_{e}^{2}}


$\displaystyle E_{\mathrm{ke}} = \frac{3}{2 m} \cdot \frac{1}{r_{e}^{2}}

An additional restraining radial potential may also be needed to limit the proton size.

Then, modify the function ElecElecPauli so that the interaction does not apply to non-classical nuclei. Finally, change functions such as ProcessInputFile(), ProcessNucleus(), and ProcessElectrons() in eff_input to read in non-classical nuclei from the input file; functions in eff_output would need to be changed too.

Create and destroy particles dynamically

To create and destroy particles dynamically, we first need to change the data structure used to hold nucleus and electron data from a fixed array to a linked list.

To do this, redeclare the variables nuc and elec in eff_global as type LIST_PTR:

  LIST_PTR *nuc, *elec

Rewrite the functions Initialize(), AddNucleus(), and AddElectron() in eff_create to work with the linked lists, using the functions new_list_ptr() and add_list_ptr() from eff_listptr. Also, use clear_list_ptr() and delete_list_ptr() to clear and deallocate the lists at the end of the program.

Next, modify the functions UpdateEnergyForces(), UpdateKineticEnergy(), UpdatePauliPeriodic(), UpdateWidthRestraint(), and AssignParticleBins() in eff_update so that they iterate over the linked list instead of over the array elements. For example, a loop to iterate over the nuclei:

  for (ptr = nuc->head; ptr != 0; ptr = ptr->next)
    NUC *nuc_ptr = (NUC *) ptr->element;
The functions in eff_ewald will need to be updated in the same way too.

Change the accessor functions in eff_access so instead of taking a nuclear or electron index as input, they take a pointer from the linked list. All the functions that reference eff_access functions will then need to be changed to loop over linked-list iterators, e.g. functions in eff_dynamics and eff_minimize. The alternative is to preserve the indexes, but to map the indexes onto linked list elements.

Finally, add a function to eff_listptr to delete elements from a linked list. Now it should be possible to create and delete particles in the middle of a simulation.

There are other issues to consider as well, such as energy conservation, making sure new particles do not overlap old ones too much, etc. Also when outputting the particle positions, energies, etc. in the .eff output file, it would be useful for trajectory analysis to be able to keep track of which particles were the same in between frames. For this purpose it would be useful to associate particles with unique indices or tags that would be printed out along with the per-particle information.

next up previous contents index
Next: To access hidden features Up: eFF Programmers' Guide Previous: To modify current features   Contents   Index
Julius 2008-04-29