Previous Thread
Next Thread
Print Thread
Page 1 of 3 1 2 3
electron density profiles via COOR HIST Z
#7329 06/28/05 08:42 PM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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

Re: electron density profiles via COOR HIST Z
rmv #7330 07/03/05 04:10 PM
Joined: Jan 2005
Posts: 86
M
Forum Member
Offline
Forum Member
M
Joined: Jan 2005
Posts: 86
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

Re: electron density profiles via COOR HIST Z
mzzt #7331 07/03/05 04:34 PM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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).


Rick Venable
computational chemist

Re: electron density profiles via COOR HIST Z
rmv #7332 07/03/05 05:19 PM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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).


Rick Venable
computational chemist

Re: electron density profiles via COOR HIST Z
rmv #7333 07/06/05 01:30 PM
Joined: Jan 2005
Posts: 86
M
Forum Member
Offline
Forum Member
M
Joined: Jan 2005
Posts: 86
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.

Re: electron density profiles via COOR HIST Z
mzzt #7334 07/06/05 01:34 PM
Joined: Jan 2005
Posts: 86
M
Forum Member
Offline
Forum Member
M
Joined: Jan 2005
Posts: 86
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.

Re: electron density profiles via COOR HIST Z
rmv #7335 05/01/08 10:32 PM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
Posts: 165
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

Re: electron density profiles via COOR HIST Z
beginner #7336 05/01/08 10:49 PM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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.

Last edited by rmv; 02/20/09 02:00 AM.
Re: electron density profiles via COOR HIST Z
rmv #7337 05/01/08 11:12 PM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
Posts: 165
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

Re: electron density profiles via COOR HIST Z
beginner #7338 05/02/08 12:05 AM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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

Re: electron density profiles via COOR HIST Z
rmv #7339 05/30/08 07:44 AM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
Posts: 165
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

Re: electron density profiles via COOR HIST Z
rmv #7340 05/30/08 04:37 PM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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

Re: electron density profiles via COOR HIST Z
rmv #7341 06/25/08 10:38 AM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
Posts: 165
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?

Re: electron density profiles via COOR HIST Z
beginner #7342 06/25/08 05:46 PM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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.

Re: electron density profiles via COOR HIST Z
rmv #7343 06/26/08 08:55 AM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
Posts: 165
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 Files
18501-lip_sod.jpg (0 Bytes, 684 downloads)
Last edited by beginner; 06/26/08 09:56 AM.
Re: electron density profiles via COOR HIST Z
beginner #7344 06/26/08 10:55 AM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
Posts: 165
the plot for water indicates some errors when i integrate from 0 to -z

Attached Files
18502-water_charg.jpg (0 Bytes, 713 downloads)
Re: electron density profiles via COOR HIST Z
beginner #7345 06/26/08 04:08 PM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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.

Re: electron density profiles via COOR HIST Z
rmv #7346 06/27/08 03:12 PM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
Posts: 165
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 Files
18524-lip_sod_short_charg_poup1.jpg (0 Bytes, 689 downloads)
Re: electron density profiles via COOR HIST Z
rmv #7347 02/09/09 01:21 PM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
Posts: 165
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

Re: electron density profiles via COOR HIST Z
beginner #7348 02/09/09 11:54 PM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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.

Re: electron density profiles via COOR HIST Z
rmv #7349 02/12/09 04:12 PM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
Posts: 165
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

Re: electron density profiles via COOR HIST Z
beginner #7350 02/12/09 11:06 PM
Joined: Sep 2003
Posts: 8,580
Likes: 11
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,580
Likes: 11
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.

Page 1 of 3 1 2 3

Moderated by  lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.4
(Release build 20200307)
Responsive Width:

PHP: 7.3.31-1~deb10u1 Page Time: 0.034s Queries: 61 (0.029s) Memory: 0.8429 MB (Peak: 0.9949 MB) Data Comp: Off Server Time: 2021-11-30 23:44:14 UTC
Valid HTML 5 and Valid CSS