Two separate scripts follow; the first computes the total electron density, the second that of the lipids only. The time interval for averaging is controlled by the number of files processed, and an additional argument for the maximum Z is expected. These scripts are run via e.g.

! SET NO. OF ELECTRONS BASED ON ELEMENT TYPE scalar wmain set 1. sele type H* end scalar wmain set 6. sele type C* end scalar wmain set 7. sele type N* end scalar wmain set 8. sele type O* end scalar wmain set 15. sele type P* end scalar wmain store 1

! SETUP TO READ FILES set kfi @FFI set nfrm 1 ! THIS IS THE LOOP OVER FRAMES label lp1

if nfrm gt 1 goto noway ! OPEN A FILE WHEN NFRM IS SET TO 1 open unit 21 read file name dyn@KFI.trj traj iread 21 nread 1 label noway

! READ A FRAME, DO AN IMAGE UPDATE traj read update imgfrq 1

! PER FRAME SLAB VOLUME; ALLOWS FOR NPgT calc vslb = ?XTLA * ?XTLB * @ZIN scalar wmain recall 1 scalar wmain divi @VSLB

! TEST FOR FRAMES PER FILE incr nfrm by 1 if nfrm le ?NFILE goto lp1 set nfrm 1 close unit 21 ! TEST FOR SEQUENTIAL FILE NO. incr kfi by 1 if kfi le @LFI goto lp1

bomlev -2 open unit 12 write card name tot@LFI.dat coor hist Z hnum @NBN hmin @ZLM hmax @MXZ sele none end hprint - hnorm ?NCONFIG iunit 12

stop

Since there's no atom selection for the COOR HIST Z command in the loop, this computes the total electron density. The next script does just the lipids--

* slab e- density along Z in intervals; max Z = @MXZ * lipids only; first file @FFI last file @LFI *

! SET NO. OF ELECTRONS BASED ON ELEMENT TYPE scalar wmain set 1. sele type H* end scalar wmain set 6. sele type C* end scalar wmain set 7. sele type N* end scalar wmain set 8. sele type O* end scalar wmain set 15. sele type P* end scalar wmain store 1

! SETUP TO READ FILES set kfi @FFI set nfrm 1 ! THIS IS THE LOOP OVER FRAMES label lp1

if nfrm gt 1 goto noway ! OPEN A FILE WHEN NFRM IS SET TO 1 open unit 21 read file name dyn@KFI.trj traj iread 21 nread 1 label noway

! READ A FRAME, DO AN IMAGE UPDATE traj read update imgfrq 1

! PER FRAME SLAB VOLUME; ALLOWS FOR NPgT calc vslb = ?XTLA * ?XTLB * @ZIN scalar wmain recall 1 scalar wmain divi @VSLB

coor hist Z hnum @NBN hmin @ZLM hmax @MXZ hsave weight - sele resn dppc end

! TEST FOR FRAMES PER FILE incr nfrm by 1 if nfrm le ?NFILE goto lp1 set nfrm 1 close unit 21 ! TEST FOR SEQUENTIAL FILE NO. incr kfi by 1 if kfi le @LFI goto lp1

bomlev -2 open unit 12 write card name edens/lpd@LFI.dat coor hist Z hnum @NBN hmin @ZLM hmax @MXZ sele none end hprint - hnorm ?NCONFIG iunit 12

stop

Changing the atom selection to SELE RESN TIP3 END would compute the electron density of the water, etc.

Dear Prof. Rick, I have used these scripts for evaluating electron density of my system, but since it is an NPAT system I have removed the following lines from the script because the script stopped working if these lines were incorporated.

! PER FRAME SLAB VOLUME; ALLOWS FOR NPgT calc vslb = ?XTLA * ?XTLB * @ZIN scalar wmain recall 1 scalar wmain divi @VSLB

It is incorrect to remove those lines w/o making other changes; they are needed for all ensembles, NVE, NVT, NPAT, or NPgT the way the script is currently written. By dividing by the slab volume for each stored frame, volume changes are allowed, and the method is not limited to constant volume or constant surface area simulations. If you remove those lines, you must divide by the slab volume at some other point (and you can't evaluate NPgT simulations).

Dear Prof Rick, Thank you so much for your help. I got the vslab lines incorporated in the script and have got the results.

Some of the points in the script are still not clear to me, e.g.,

1. The "scalar keyname set" value gives some numbers to the different selections i.e. scalar wmain set 1. sele type H* end scalar wmain set 6. sele type C* end scalar wmain set 7. sele type N* end scalar wmain set 8. sele type O* end scalar wmain set 15. sele type P* end

The numbers are in the increasing order. Is it according to the charge on the atoms since the question here is to calculate the electron density? So, P being negative gets the highest number and O,N and C ~ similar and H the least i.e a value of 1?

2. Does the following line

"scalar wmain store 1"

stores the numbers (and hence the electron density) in wmain and stores this for H, C, N, O and P in the store number 1?

Eagerly waiting for your reply, Regards, Priyanka.

Dear Prof Rick, I am extremely sorry for my previous post. It just occured to me after sending my previous post that these "set" numbers corresponds to the number of electron in the atoms selected. Sorry to bother you Sir, Apologies, Priyanka.

sorry if this is a naive question. I am going to modify this script to calculate electrostatic potential for the membrane. However, I am not sure how to make a double integration based on your script.

The script can be used to get the charge density, which can be saved to a file. Just put the charges into WMAIN instead of the atomic number; see scalar.doc

The double integral is a separate step; a simple program in the language of your choice may be needed for that. I chose Fortran.

Unless you can modify and compile Fortran source code, my program may not be of much use. If you plan to work in a comp chem field, you should learn a programming language. I've heard many good things about the Python language, which has numeric and biomolecule related packages available.

There are commands in CORREL, which might be okay for this; the charge density file can be read into a time series array via READ, and the time series manipulation commands include an integration command. Try something like--

! ASSUMES 72 A Z RANGE, BINNED TO 0.1 A INTERVAL correl maxtim 720 noupdate ! DECLARE SERIES OF TYPE ZERO enter chg zero ! OPEN AND READ CHARGE FILE; SKIP*DELTA SHOULD BE THE BIN WIDTH ! OFFSET IS THE Z COORD OF THE FIRST DATA POINT open unit 1 read card name zchg.dat read chg unit 1 dumb colu 2 skip 1 delta 0.1 offset -36.0 ! INTEGRATE TWICE mantim chg integ mantim chg integ ! WRITE RESULT open unit 2 write card name chgint2.dat write chg unit 2 dumb time * dummy * end ! EXIT CORREL

I've had some trouble with using an 'edit-spec' with that READ command myself lately. An alternate approach for the first part is--

! DECLARE SERIES OF TYPE ZERO enter chg zero ! SETUP TO MATCH Z DATA; SKIP*DELTA SHOULD BE THE BIN WIDTH ! OFFSET IS THE Z COORD OF THE FIRST DATA POINT edit chg skip 1 delta 0.1 offset -36.0 ! OPEN AND READ CHARGE FILE open unit 1 read card name zchg.dat read chg unit 1 dumb colu 2

thank you very much. I obtain an reasonable graph when I did integrate from z to +z. However, the result seemed wrong if I integrate from -z to z. do you have any suggestion to modify the input or only thing i can do is reorganize my data?

The result depends on whether the integration starts from the middle of the bilayer or the middle of the water; the sign changes. The starting region for the integration should also be fairly flat, or else the integral will have an apparent slope.

Opps! I did a mistake, I wanted to say that I got reasonable graph when I integrate from 0 to +z, but not from -z to 0. I would expect a symmetric plot (from -z to 0 and from 0 to +z). And I also expect a fairly flat at the stating region before the slope appears, but i didnt get it (although the charge density show zero value at that region). Please find the enclosed picture. Any comment on that? thank you very much

For the slope problem, try shifting the starting point for the integration, perhaps by dropping a few initial points. I don't have much idea about the problem with the water. I've generally done the double integration using the -z to +z range.

Also, you need to average the charge density over many ns of data.

wrt to the script that we discuss, it seemed to me that I need to make an integration separately, from 0 to +z and from 0 to -z. if I did the double integration from -z to +z , I got a strength output as you can see from the attachment.

CHARMM uses a distinct system of units, the AKMA system. I.e. Angstroms, Kilocalories/Mole, Atomic mass units. All distances are measured in Angstroms, energies in kcal/mole, mass in atomic mass units, and charge is in units of electron charge.

but in order to get the electrostatic potential one said by integrating the poisson equation over the box in the z coordinate : Integral(0-z)dz´ Integral(0-z´)[ Q(z´´) dz´´ ]. I am wondering if the INTEGRATE TWICE as you suggested give the same meaning as electrostatic potential.

In addition, this parameter would be depend on the force field. I am right?

The double integral of the charge density across the membrane gives the dipole potential; I'm not quite sure if that's what you want or not. It is strongly dependent on the force field.

By default CHARMM assumes a time base, so the stated definition of INTEgrate is wrt. to time starting at zero. For a z profile, you have to change the independent variable to match the z values via EDIT; for a range of -40 to 40 in 0.1 A intervals, something like

edit chg delta 0.1 skip 1 offset -40.

would be needed.

Perhaps you should reconsider writing your own program for this.