CHARMM Development Project
Posted By: rmv electron density profiles via COOR HIST Z - 06/28/05 08:42 PM
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.
Posted By: mzzt Re: electron density profiles via COOR HIST Z - 07/03/05 04:10 PM
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
Posted By: rmv Re: electron density profiles via COOR HIST Z - 07/03/05 04:34 PM
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).
Posted By: rmv Re: electron density profiles via COOR HIST Z - 07/03/05 05:19 PM
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).
Posted By: mzzt Re: electron density profiles via COOR HIST Z - 07/06/05 01:30 PM
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.
Posted By: mzzt Re: electron density profiles via COOR HIST Z - 07/06/05 01:34 PM
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
Posted By: rmv Re: electron density profiles via COOR HIST Z - 05/01/08 10:49 PM
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
Posted By: rmv Re: electron density profiles via COOR HIST Z - 05/02/08 12:05 AM
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
sorry, but all I got is zero
CORREL> read chg unit 12 dumb colum 2 skip 1 delta 0.2 offset -34.9
**** Warning **** The following extraneous characters were found while command processing in CORREL
SKIP 1 DELTA 0.2 OFFSET -34.9

NSER: 1
NAMES: CHG
TOTALS: 350
AVERAGE: 0.303590
FLUCT: 0.098488
VECCOD: 1
CLASS: ZERO
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 1
VALUE: 0.000000

CORREL> mantime chg integ
NSER: 1
NAMES: CHG
TOTALS: 350
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: ZERO
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 1
VALUE: 0.000000

CORREL> mantime chg integ
NSER: 1
NAMES: CHG
TOTALS: 350
AVERAGE: 0.000000
FLUCT: 0.000000
VECCOD: 1
CLASS: ZERO
VELCOD: 0
SKIP: 1
DELTA: 0.000000
OFFST: 0.000000
GECOD: 1
VALUE: 0.000000

any comments on that?
thank you very much
Posted By: rmv Re: electron density profiles via COOR HIST Z - 05/30/08 04:37 PM
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
dear rick,

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?
Posted By: rmv Re: electron density profiles via COOR HIST Z - 06/25/08 05:46 PM
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

Attached picture 18501-lip_sod.jpg
the plot for water indicates some errors when i integrate from 0 to -z

Attached picture 18502-water_charg.jpg
Posted By: rmv Re: electron density profiles via COOR HIST Z - 06/26/08 04:08 PM
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.

Attached picture 18524-lip_sod_short_charg_poup1.jpg
Quote:

The script can be used to get the charge density, which can be saved to a file.






Dear rick, what is the unit of the y-axis in this case? is it coulomb/cubic angstrom or volt?

thank you
Posted By: rmv Re: electron density profiles via COOR HIST Z - 02/09/09 11:54 PM
From usage.doc--

The CHARMM system of units: AKMA.

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.
Quote:


! INTEGRATE TWICE
mantim chg integ
mantim chg integ
! WRITE RESULT





form the correl.doc

INTEgrate Q(z) = Integral(0-z) [ Q(z) dz ]

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?

thank you
Posted By: rmv Re: electron density profiles via COOR HIST Z - 02/12/09 11:06 PM
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.
© CHARMM forums