Calculate electrostatic potential of a bilayer with respect to z This wiki explains how to calculate the total electrostatic potential of a bilayer with respect to position in z. The dipole potential drop can then be found by subtracting the electrostatic potential of bulk water from the potential at the center of the bilayer. (Using the scripts provided, the electrostatic potential is defined as zero at the center of the bilayer.)
- Location of files on ZT1: /afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/epp_bilayer
- Matlab is available for free download through TERPware: https://terpware.umd.edu/Windows
1. Submit “chrg_dens.inp” using the submission file “chrg_dens.csh”. The input file “chrg_dens.inp” computes the total charge density along z (the bilayer normal) in slabs of thickness 0.1 Angstroms, from -32.95 to 32.95 Angstroms. It reads in the re-centered trajectories from the current directory with file names of the form “dyn@NFIL-recenter.dcd”. The script is set up to read inputs and trajectories from the current directory.
First update the inputs to suit your system:
- Topology file: toppar.str
- Protein structure file: step5_assembly.xplor.psf
- Coordinate file: step5_assembly.crd
- Thickness of box at start of simulation (“A”): xa
- First trajectory file number: nfil
- Last trajectory file number: In line 80, you can change “100”: “if nfil .le. 100 goto lp1”
And, in the file <chrg_dens.csh>, update the path of the CHARMM executable to your preferred CHARMM executable.
Data file: “chrg_dens.dat”
The data file produced by “chrg_dens.inp” will be called “chrg_dens.dat”. This file has two columns. The first column is the beginning of the z slab in Angstroms (of width 0.1 Angstrom) and the second column is the electron density computed in a frame. This file will be rather large.
2. Use Matlab script “epp_bilayer.m” to process data file “chrg_dens.dat”.
The script “epp_bilayer.m” symmetrizes the electron density, integrates the electron density with respect to z, converts to volts, establishes electrostatic potential = 0 at the bilayer center, and plots the electrostatic potential as a function of z. It was written to process a pure DPPC bilayer, but will work with any system that is symmetric along z.
The vector containing the electrostatic potential relative to the bilayer center is “epp_dppc1“, and units of this vector are Volts.