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

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.011s Queries: 34 (0.006s) Memory: 0.7795 MB (Peak: 0.8593 MB) Data Comp: Off Server Time: 2021-12-01 00:11:16 UTC
Valid HTML 5 and Valid CSS