CHARMM Development Project Forums CHARMM Community Script Archive electron density profiles via COOR HIST Z

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.

charmm FFI:1 LFI:20 MXZ:40. < tot-eden.inp >& tot-eden20.out

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

! SET Z INTERVAL (SLAB THICKNESS)
set zin 0.2

calc nfi = 1 + ( @LFI - @FFI )
calc zlm = -1. * @MXZ
calc nbn = int( ( 2 *  @MXZ ) / @ZIN )

wrnlev -4

! READ RTF, PARAM, PSF, CRD FILES
stream rtfprm.str
stream psfcrd.str
stream cryst.str

! 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

! 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 Z INTERVAL (SLAB THICKNESS)
set zin 0.2

calc nfi = 1 + ( @LFI - @FFI )
calc zlm = -1. * @MXZ
calc nbn = int( ( 2 *  @MXZ ) / @ZIN )

wrnlev -4

! READ RTF, PARAM, PSF, CRD FILES
stream rtfprm.str
stream psfcrd.str
stream cryst.str

! 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

Is this right?

Eagerly waiting for your reply,
Regards,
Priyanka

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).

On further review, it is absolutely incorrect to remove ALL of those lines from that location; in particular

scalar wmain recall 1

is essential, as that command places the number of electrons for each element type in the weighting array (WMAIN).

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.

Dear rick,

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.

Thank you in advance

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.

thank you,

sorry but please ignore it if I disturb you too much.
I am not good in programing..it would be very kind of you if you can share it.

thank you

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