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.


Rick Venable
computational chemist