Page 1 of 2 1 2 >
Topic Options
#8389 - 09/29/05 11:04 PM ion addition by selective water replacement
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8389
Loc: 39 03 48 N, 77 06 54 W
The following scripts demonstrate a sampling procedure that starts with a well-minimized, solvated system. The main script, ion-emin.inp, samples 100 ion configs and keeps only the minimum energy config. The system I was using was a PC bilayer with an added anionic lipid, but the method is generally applicable to most solvated biomolecule systems. As for how many ions to add, computing concentration from the unit cell volume can be tricky, as a large portion of it may be non-water, i.e. the interior of a protein or a lipid membrane. We prefer to use the fact that water itself is ca. 55 M, and hence replacing one water of each 55 gives a 1 M solution. Adding one Na+Cl- pair for every 300 water molecules gives an ion concentration of a bit under 200 mmol, which is more like biological ionic strength. In this example, the anionic lipid has a -4 charge, so four extra Na+ ions are added to get a system with a net zero charge, which is best for Ewald methods. The script uses command line args to pass a seed for the random number generator, and to assign a run number for use in identifying output files. On a single processor, it could be run via

charmm RUN:1 SD:344568 < ion-emin.inp > ion-emin1.out

For parallel usage, which may be desirable, the script can be renamed to ion-emin.str and executed as a stream file via e.g.

* cover script for parallel; cannot rewind or backspace std input
* CHARMM loops are okay in stream files, just not stdin
*
stream ion-emin.str
stop

Other details should (hopefully) be clear from the following scripts; in addition to the main control script, the addions.str utility script is included. By setting some '@' parameters (variables) in the main script, I've used this unchanged as a utility routine in numerous projects. It's an illustration of my "modular" CHARMM programming style, which features re-usable utility scripts such as this.
_________________________
Rick Venable
computational chemist


Top
#8390 - 09/29/05 11:17 PM ion addition by selective water replacement; ion-emin.str [Re: rmv]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8389
Loc: 39 03 48 N, 77 06 54 W

* DPPC bilayer; ion placement sampling; seed @SD run @RUN
*
bomlev -1

! READ RTF AND PARAM FILES
stream rtfprm.str

! SETUP STUFF; PROJECT DEPENDENT
set icl 12 ! no. chlorides
set ina 16 ! no. sodiums
set mnd 5.5 ! minimum distance to solute, other ions
! SOLUTE EXCLUSION; selection w/o SELE ... END; SEE addions.str
set sol segid L .or. segid P
set emin 0. ! initial min energy value
set ncfg 1 ! initialize loop counter
set last 100 ! no. of passes thru the loop
! CHANGING ISEED WILL SAMPLE 100 DIFFERENT PLACEMENTS
random uniform iseed @SD
! KEEP A RECORD OF THE ENERGIES, FRAMES
open unit 9 write card name ionmc@RUN.dat

! BEGUN THE LOOP HAS
label placion
open unit 10 write card name loop@RUN.log
outu 10
! LOOP LOG FILE OVERWRITTEN ON EACH PASS; GREATLY REDUCES OUTPUT
! LAST PASS ALWAYS SAVED (IN EVENT OF FAILURE)

! FULLY SOLVATED, MINIMIZED WITH CRYSTAL LATTICE
! READ MATCHING PSF AND COOR FILES
open unit 3 read card name f45d02.psf
read psf card unit 3
close unit 3
open unit 3 read card name f45d02.crd
read coor card unit 3
close unit 3

! RANDOM WATER REPLACEMENT
stream addions.str

! ASSUMES WATER IS SEGID WAT; RENUMBER THE WATER MOLECULES
join wat renum

! SETUP CRYSTAL (DEFINE, BUILD), IMAGE CENTERING W. MODIFIED PSF
stream cryst.str

! SETUP SHAKE, NONBOND
shake bonh param
update inbfrq 5 atom vatom cutnb 12.0 ctofnb 10. cdie eps 1. -
ctonnb 8. vswitch cutim 12.0 imgfrq 5 wmin 1.0 -
ewald pmew fftx 48 ffty 48 fftz 64 kappa .33 spline order 6

! BRIEF MIN OF IONS INSERTED INTO SOLVATED MODEL
mini sd nstep 10 nprint 1
mini abnr nstep 25 nprint 5

! LOG SOME DATA
write title unit 9
* @NCFG ?ENER @EMIN
*

! KEY LOGICAL TEST; CONFIGS WITH HIGHER ENERGY REJECTED
if emin .lt. ?ENER goto test

! WRITE THE LATEST MINIMUM ENERGY RESULT; CONFIG ACCEPTED
open unit 11 write card name f45d02i@RUN.psf
write psf card unit 11
* 71 DPPC bilayer & PIP2 (-4) and 16*Na+ 12*Cl- for Ewald neutral
*

! DO AN IMAGE UPDATE AND WRITE THE COOR FILE
update
open unit 11 write card name f45d02i@RUN.crd
write coor card unit 11
* DPPC bilayer with PIP2 (-4) and 16*Na+ 12*Cl- for Ewald neutral
* run @RUN seed @SD MC @NCFG C ?XTLC E ?ENER PREV @EMIN
*

! UPDATE MINIMUM ENERGY
set emin ?ENER

! TEST FOR EXIT, AND SETUP FOR NEXT PASS; REVERT TO STDOUT, CLOSE LOG
label test
crystal free
shake off
incr ncfg by 1
outu 6
close unit 10
if ncfg le @LAST goto placion
close unit 9

return

Top
#8391 - 09/29/05 11:20 PM Re: ion addition by selective water replacement; addions.str [Re: rmv]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8389
Loc: 39 03 48 N, 77 06 54 W

* sodium and chloride addition; @ICL CL- @INA NA+ MIN DIST @MND
* excluding @SOL
*

! sol particles to exclude; atom selection w/o SELE and END
! icl no. of chlorides
! ina no. of sodiums
! mnd min distance to another solute particle

! INITIALIZE CHLORIDES; COORDS SET 0,0,0
read sequ cla @ICL
generate cl noang nodihe first none last none
coor set xdir 0.0 ydir 0.0 zdir 0.0 sele segid cl end

! PICK A WATER MOLECULE AND REPLACE IT WITH A CHLORIDE ION
set k 1
label ilp
define xcld sele ( @SOL ) .or. segid CL end ! EXCLUDED SEGMENTS
define prox sele xcld .around. @MND end ! NEARBY EXCLUDED ATOMS
define list sele atom WAT * OH2 .and. .not. prox end ! WATERS NOT NEARBY
calc in int( ?RAND * ?NSEL ) ! RANDOM INTEGER, BASED ON WATER COUNT
if in .lt. 1 set in ?NSEL ! CHECK FOR ZERO, CHANGE TO MAX VALUE
define targ sele list .subset. @IN end ! PICK WATER VIA RANDOM INDEX
coor stat sele targ end ! GET OH2 ATOM COORDS, ASSIGN TO CHLORIDE
coor set xdir ?XAVE ydir ?YAVE zdir ?ZAVE sele atom CL @K CLA end
delete atom sele .byres. targ end sort ! REMOVE THE WATER MOLECULE
incr k by 1
if k .le. @ICL goto ilp

! INITIALIZE SODIUMS, REPEAT ABOVE PROCESS
read sequ sod @INA
generate na noang nodihe first none last none
coor set xdir 0.0 ydir 0.0 zdir 0.0 sele segid na end

! NOTE THAT BOTH SEGID CL AND SEGID NA ARE 'EXCLUDED'
! NEARBY WATERS ARE REALLY EXCLUDED; REPLACE ONLY 'BULK' WATERS
set k 1
label klp
define xcld sele ( @SOL ) .or. segid CL .or. segid NA end
define prox sele xcld .around. @MND end
define list sele atom WAT * OH2 .and. .not. prox end
calc in int( ?RAND * ?NSEL )
if in .lt. 1 set in ?NSEL
define targ sele list .subset. @IN end
coor stat sele targ end
coor set xdir ?XAVE ydir ?YAVE zdir ?ZAVE sele atom NA @K SOD end
delete atom sele .byres. targ end sort
incr k by 1
if k .le. @INA goto klp

return

Top
#8392 - 02/24/06 11:11 AM Re: ion addition by selective water replacement; addions.str [Re: rmv]
cneale Offline
Forum Member

Registered: 10/06/05
Posts: 11
Incredibly useful script.
Thank you for the excellent post.

Top
#8393 - 04/18/08 05:37 AM Re: ion addition by selective water replacement; addions.str [Re: cneale]
zee Offline
Forum Member

Registered: 07/05/04
Posts: 53
Thank your for the very clear and well commented script.

Best regards,
zee

Top
#23857 - 04/02/10 02:29 AM Re: ion addition by selective water replacement; addions.str [Re: zee]
Kellen Offline
Forum Member

Registered: 01/02/10
Posts: 15
Loc: Tampa, FL
This script worked very well and I was able to get a system with ions that are nicely spaced apart from each other (and my protein), but I was wondering, how may I adapt this to make sure that the ions aren't placed too close to the edges of my water box? I am running MD using periodic boundary conditions (image/crystal) and noticed that one of my ions hopped to the other side of the system and is now interacting too closely with an ion that was placed close to that edge of the box.

Thanks smile
_________________________
- Kellen -

Top
#23863 - 04/02/10 02:33 PM Re: ion addition by selective water replacement; addions.str [Re: Kellen]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8389
Loc: 39 03 48 N, 77 06 54 W
The ions are fairly mobile, and will adapt to the field of the less mobile charges within a few 100 ps of simulation.

The ion placement loops would have to be modified to exclude some more water molecules, based on their coordinate values.

The energy evaluation in the calling script uses PBC, so if the placement caused a major energetic problem, the config would not have been the lowest energy coordinate set chosen.

If you think the starting config may be a problem, one option is to change the random seed and run the ion addition again. I often do 5-10 runs with different random seeds, and compare the final coordinate sets produced. If you have the computer power available, it's also not a bad idea to run simulations in triplicate, each with a different initial atom packing.
_________________________
Rick Venable
computational chemist


Top
#25419 - 09/19/10 03:01 AM Re: ion addition by selective water replacement; addions.str [Re: rmv]
Todd Offline
Forum Member

Registered: 07/12/08
Posts: 12
Thanks for sharing it.
It's awesome

Top
#27241 - 04/18/11 02:05 PM Re: ion addition by selective water replacement; addions.str [Re: Todd]
indian_user Offline
Forum Member

Registered: 03/04/11
Posts: 42
Hi,

could you attach an example of the file 'rtfprm.str' ?

thanks.

Top
#27246 - 04/18/11 03:01 PM Re: ion addition by selective water replacement; addions.str [Re: indian_user]
lennart Offline

Forum Member

Registered: 09/25/03
Posts: 4743
Loc: ~ 59N, 15E
possibly somehting like
* my little file that reads my RTF and PARAmeter files
*
read rtf card name my_favorite_rtf.rtf
read para card name my_favorite_param.inp
return
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
Page 1 of 2 1 2 >

Moderator:  chmgr, John Legato, petrella