CHARMM Element doc/pimplem.doc $Revision: 1.1 $#File: PIMPLEM, Node: Top, Up: (chmdoc/perturb.doc), Next: Description Implementation of the Thermodynamic Simulation Method * Menu: * Description:: How Chemical Perturbation works.* File Formats:: Output File Formats for Chemical Perturbation.* IC Implementation:: Implementation and File Formats for Internal Coordinate Perturbation#File: PIMPLEM, Node: Description, Up: Top, Next: File Formats, Previous: Top How the Chemical Perturbation Energy Calculation Works For thermodynamic perturbation calculations the atoms making upthe system described by the hybrid Hamiltonian, H(lambda), can be dividedinto four groups. 1) The environment part - all atoms that do not changeduring the perturbation. E.g., for ethanol -> propane the solvent andthe terminal methyl group. 2) The reactant atoms - the atoms that arepresent at lambda = 0 and absent at lambda = 1. 3) The product atoms -the atoms that are absent at lambda = 0 and present at lambda = 1. 4)The COLO atoms - atoms that are present in both the reactant and productbut change charge in going from one to the other. Certain basic premises underly our approach. Energy values arefactored by lambda (or functions thereof), never the energy functionsthemselves. The standard energy routines are called unchanged and can bemodified without requiring changes to the perturbation routines as longas the calling sequence remains the same. Potential energy terms arewritten to output during a trajectory and in the case of the windowmethod trajectories can be combined. Futhermore any lambda -> lambda'can be calculated post priori and additional lambda points can be addedas desired. Most other implementations do not appear to allow this.There is, however, a price entailed namely a certain amount of redundantcalculation. Furthermore , purely as a matter of conceptual preference,the entire perturbation part of the Hamiltonian is facter by lambda inthe same way. There has been some advocacy of factoring the attractiveand repulsive part of the Lennard-Jones potential with different powersof lambda (see Cross). We want to calculate the potential energy U(lambda) = Uenv +(1-lambda)**N Ureac + lambda**N Uprod, where N is positive integerexponent and Uenv is the energy of the common environment part of thesystem. The residue topology file for the system undergoing theperturbation has all the internal coordinate terms for both the productand reactant parts and the regular CHARMM energy routine calculates anenergy term that in it's sum contains part of Ureac and Uprod along withUenv and in certain cases, as will be discussed shortly, an additionalterm that needs to be removed. The residue description must containnon-bonded exclusions between the product and reactant atoms. Of course,none of this is factored correctly, or at all, by lambda. The approach to obtaining a the correct U(lambda) is an indirectone. Instead of making it so that the normal energy routine calculatesUenv only and having the perturbation energy routine calcuated determine(1-lambda)**N Ureac + lambda**N Uprod, we have it instead calculate theamount that must be subtracted from the normal energy routine value (hereafter referred to as Unorm) to get U(lambda). The previous statementmust be amended for the case where there are COLO atoms. Then, Unormcontains a term that must be totally removed and is missing some termscompletely, which must be added. For the internal coordinate energy terms and the non-bonded vander Waals interactions, the amount that must be subtracted from Unorm toobtain U(lambda) is given by: U(lambda) = Unorm + Ucorr since, U(lambda) = U(env) + (1 - lambda)**N Ureac + lambda**N Uprodand Unorm = U(env) + Ureac + Uprodthen -Ucorr = [1-(1-lambda)**N]Ureac + [1-lambda**N]Uprod . We have currently ignored the electrostatic terms. If there are no COLOatoms the above expressions hold true for those terms as well. If there are COLO atoms , the situation becomes a bit morecomplicated. To discuss this the following nomenclature is introduced: [reac| reac,colo-r,env] The expression above indicates the calculation of the electrostaticinteraction between reactant atoms and 1) other reactant atoms 2) COLOatoms with the reactant energy charges and 3) with environment atoms. Unorm contains the following electrostatic terms: [reac| reac, colo-r, env] + [color | prod, colo-r, env] + [prod | prod, env] The term [ colo-r | prod ] must be removed in it's entirety (productatoms do not interact with reactant charges (colo-r). And the missinginteractions involving colo-p (product) charges must be added (suitablyfactored by lambda). To do this Ucorr must contain: (1 - (1 - lambda)**N) { [reac | reac, colo-r, env] + [colo-r | colo-r, env] } + (1 - lambda**N)[prod | prod, env] + 1[color|prod] - lambda**N [colo-p | colo-p, prod, env] Note that -Ucorr is passed from the perturbation energy routine, thus thenegative term (last one) actually adds what is totally missing fromUnorm. The electrostatic contribution to Ucorr is actually calculated inan even more round-about fashion than that which is given above. First both the van der Waal's and electrostatic interactionsinvolving reactant and product atoms with everything (except interactionsbetween reactant and product atoms) are calculated. The reactant colo-rcharges are used for this. This provides the term: (1 -(1-lambda)**N)[reac | env, colo-r, reac ]and (1 - lambda**N)[prod | env, colo-r, prod ]. If there are no COLO atoms, this is all we need (absent the colo-r termin the expressions). Otherwise, three more calculations, involving onlythe electrostatic energy, are required. The first involves interactionsbetween colo-r charges with environment and other colo-r charges: (1-(1-lambda)**N)[colo-r | env, colo-r] Next the colo-r product atom electrostatic interaction is calculated andfactored by a function of lambda that compensates for the amount in thesecond ([prod | colo-r ...] ) calculation. In that term,1-lambda**N[prod | colo-r] is included so we must determine, (lambda**N)[colo-r| prod]