Questions and problems?

LEaP

        Why does the solvatebox command result in an initial density
        of ~0.74 g/cc?

Superimposing a box of water on the solute & box and then subtracting overlapping waters leads to vdw voids at the solute & box interfaces, since perfect packing is not achieved. Thus one should always choose a larger box size than desired when using the default closeness parameter of 1.0, and use constant pressure at some point during equilibration to achieve proper density. As a rule of thumb, add 2.5A to the desired boxsize. A denser initial system can be achieved using closeness < 1.0; in this case there will be relatively close contacts of solute and solvent.

> Whenever I load a .pdb file into xLEaP, it generates numerous new atoms
> However, they are generated without a type, charge, pert.name, and
> pert.type. Could someone please help me in regard to what can be done to
> prevent that, or what to do when those 'extra' atoms are generated.

When LEaP loads a PDB file it tries to match the residue name in the PDB file to a residue in the LEaP database, then tries to match all the atoms in the PDB residue to the LEaP prototype residue.

For example, a SER residue:

ATOM      1  N   SER     0       3.326   1.548  -0.000  1.00  0.00
ATOM      2  H   SER     0       3.909   0.724  -0.000  1.00  0.00
ATOM      3  CA  SER     0       3.970   2.846  -0.000  1.00  0.00
ATOM      4  HA  SER     0       3.672   3.400  -0.890  1.00  0.00
ATOM      5  CB  SER     0       3.577   3.654   1.232  1.00  0.00
ATOM      6  HB2 SER     0       2.497   3.801   1.241  1.00  0.00
ATOM      7  HB3 SER     0       3.877   3.116   2.131  1.00  0.00
ATOM      8  OG  SER     0       4.231   4.925   1.197  1.00  0.00
ATOM      9  HG  SER     0       3.983   5.434   1.973  1.00  0.00
ATOM     10  C   SER     0       5.486   2.705  -0.000  1.00  0.00
ATOM     11  O   SER     0       6.009   1.593  -0.000  1.00  0.00

LEaP reads this PDB residue in, sees that it has a SER residue, and all the atom names match perfectly to the LEaP prototype.
If a residue in the pdb file has atoms with names that do not exactly match the atom names that leap knows, it will assume that those atoms are missing and create new ones so that the residue has all its atoms. The atoms in the pdb will not match atoms in leap's databases and their types and charges will be unknown.
If the pdb file has a residue that leap does not find in its database, it will not be able to assign atom tyoes or charges, etc since they simply are unknown to leap. If I were to add a new atom with a new name to the SER residue in the PDB file, that atom would not appear in the LEaP prototype residue, and leap would get confused.

There are several solutions:

"Massage" the input PDB file to have the "correct" atom and residue names, with appropriate TER records between molecules.
Use a PDB atom map, in LEaP (see "help addpdbatommap" and the default leaprc)

which residues give extra atoms when loading the original PDB file, then resolve the problems using a copy of the original PDB file.

LEaP can't find some new atom types I added for my molecule. BUT both in the EDIT the unit and EDIT SELECTED ATOMS both atom OY & HY are present:

	> saveAmberParmPert solvate sol.top sol.crd 
	Checking Unit.  
	Building topology.  
	Building atom parameters.  
	For atom: .R[cpa 1].A[H3 8] Could not find type: HY 
	For atom: .R[cpa 1].A[H3 8] type OY
		Could not find perturbed type:  HY
	Parameter file was not saved.

Atoms have names and types. Names are unique within the residue, and serve to distinguish the individual atoms (their topological positions) - e.g. H1'. Types are not unique - they are like a generalization of element to include (encode) information such as hybridization and bonded neighbors - e.g. amino hydrogen. The force field is essentially a definition of the properties of different atom types. These are all defined in the parm94.dat (new ff) or parm91X.dat (old ff). If you add a new type, you must define it - this is normally done in a frcmod file. However, one does not want to add new atom types to the force field unless absolutely necessary, since it means also generating bond, angle, and torsion parameters for the new type. I suspect that, especially in your case of hydrogen (HY), there may already an appropriate type that you could choose from the force field for the atom in its perturbed state.

Please explain the output from 'solvatebox':


	>solvateBox A WATBOX216{2 0.1 0.1 }

	  Solute vdw bounding box:              12.686 9.519 7.933
	  Total bounding box for atom centers:  16.686 9.719 8.133
	  Solvent unit box:                     18.774 18.774 18.774
	  Total vdw box size:                   18.743 11.741 11.491 angstroms.

Explanation:


	  Solute vdw bounding box:              12.686 9.519 7.933

This box contains the solute atoms - including vdw.


	  Total bounding box for atom centers:  16.686 9.719 8.133

This contains the solute atom centers plus twice the distances {2 0.1 0.1 } since solvent is placed _around_ the solute.


	  Solvent unit box:                     18.774 18.774 18.774

Size of the box of 216 water molecules (WATBOX216) that will be used to overlay the solute, remove water molecules that are too close to or inside the solute, and be trimmed to the desired size.


	  Total vdw box size:                   18.743 11.741 11.491 angstroms.

Total bounding box for solute+solvent atom centers plus 2 * radius of water. This is the 'final' box that sander/gibbs/roar will see in prmtop.



	  Total mass 341.320 amu,  Density 0.224 g/cc
	  Added 13 residues.

I am defining several atom types unique to retinoic acid... and ... I have prepared a parameter file that contains the masses of my atom types along with BOND, ANGL, DIHE, and NONB parameters. ...

When I try to read this parameter file into Leap (using the loadamberparams command), I get an error message, a part of which is shown below:


Unknown keyword: HCR -CR10- in parameter file.
Unknown keyword: HCR -CR11- in parameter file.
Unknown keyword: HCR -CR11- in parameter file.
Unknown keyword: HCR -CR12- in parameter file.

Atom types in Amber are at most two characters long, not the three or four characters in the example above. Note that an angle definition, for example, has to be:


xx-yy-zz     force-const   ideal-angle

i.e. the beginning of the line has two characters for each atom type, with a hyphen in between, and no "extra" characters (or spaces). See, e.g. p. 122 of Part II of the manual for detailed angle formats.

Since Leap couldn't understand your angle parameters, it assumed it was seeing a keyword line that it did not understand....

I need to use the "bond" rather than "bondbydistance" command. The protocol is : bond for a single bond between atom1 and atom2. However, if I type bond C56 C57 I'm told bond: Argument #1 is type String must be of type: [atom]

How exactly is [atom] designated if not by atom name?

If the unit is X and the residue is the 1st residue in the unit (e.g. the only residue), it would be


> bond X.1.C56 X.1.C57

And in general this is how atoms are defined in leap commands.

When I asked for saveamberparmpert .... the result was error messages because of MISSING PROPER TORSION PARAMETERS:

    *** Proper torsion parameters missing ***
     atom names: N1-C7-C8-H8
     atom types: NZ-CX-CX-HZ =pert=> OZ-DC-DC-DH 
    Please add a dummy parameter of multiplicity 3 
    for the pert types to your parameter set. 
     - e.g. OZ-DC-DC-DH 1 0.0 0.  3.  
    (This is because multiple torsional potentials may apply to a
    single torsion, and each is perturbed individually in gibbs.)

When leap is complaining about missing proper torsion parameters, the thing to do is _add_ parameters of the given multiplicity, normally the ones suggested (i.e. they will contribute 0 energy). Then there might be a second round of such msgs when you try to saveamberparmpert again. In this case, you _add_ some more parameters as before. Now saveamberparmpert should work. Do not _change_ any parameters in response to these msgs unless you understand why the params are there.

I was trying to find out how to define the PMF-type perturbation option in LEaP, but I could not find it.

I have a dipeptide whose peptide bond will be changed from 0 to 180 degrees during a GIBBS run.

You should be able to

> saveamberparmpert molecule nullpert.top nullpert.crd

_without_ first defining any per-atom perturbations, and then define your PMF perturbation in the gibin file. This corresponds with the instruction for setting up a null perturbation topology file using parm, given in gibbs.doc's PMF section.