The following scripts demonstrate a sampling procedure that starts with a well-minimized, solvated system. The main script, ion-emin.inp, samples 100 ion configs and keeps only the minimum energy config. The system I was using was a PC bilayer with an added anionic lipid, but the method is generally applicable to most solvated biomolecule systems. As for how many ions to add, computing concentration from the unit cell volume can be tricky, as a large portion of it may be non-water, i.e. the interior of a protein or a lipid membrane. We prefer to use the fact that water itself is ca. 55 M, and hence replacing one water of each 55 gives a 1 M solution. Adding one Na+Cl- pair for every 300 water molecules gives an ion concentration of a bit under 200 mmol, which is more like biological ionic strength. In this example, the anionic lipid has a -4 charge, so four extra Na+ ions are added to get a system with a net zero charge, which is best for Ewald methods. The script uses command line args to pass a seed for the random number generator, and to assign a run number for use in identifying output files. On a single processor, it could be run via
charmm RUN:1 SD:344568 < ion-emin.inp > ion-emin1.out
For parallel usage, which may be desirable, the script can be renamed to ion-emin.str and executed as a stream file via e.g.
* cover script for parallel; cannot rewind or backspace std input
* CHARMM loops are okay in stream files, just not stdin
Other details should (hopefully) be clear from the following scripts; in addition to the main control script, the addions.str utility script is included. By setting some '@' parameters (variables) in the main script, I've used this unchanged as a utility routine in numerous projects. It's an illustration of my "modular" CHARMM programming style, which features re-usable utility scripts such as this.