Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
ion addition by selective water replacement
#8389 09/30/05 03:04 AM
Joined: Sep 2003
Posts: 8,503
rmv Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 8,503
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

ion addition by selective water replacement; ion-emin.str
rmv #8390 09/30/05 03:17 AM
Joined: Sep 2003
Posts: 8,503
rmv Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 8,503

* 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

Re: ion addition by selective water replacement; addions.str
rmv #8391 09/30/05 03:20 AM
Joined: Sep 2003
Posts: 8,503
rmv Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 8,503

* 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

Re: ion addition by selective water replacement; addions.str
rmv #8392 02/24/06 04:11 PM
Joined: Oct 2005
Posts: 11
C
Forum Member
Offline
Forum Member
C
Joined: Oct 2005
Posts: 11
Incredibly useful script.
Thank you for the excellent post.

Re: ion addition by selective water replacement; addions.str
cneale #8393 04/18/08 09:37 AM
Joined: Jul 2004
Posts: 53
zee Offline
Forum Member
Offline
Forum Member
Joined: Jul 2004
Posts: 53
Thank your for the very clear and well commented script.

Best regards,
zee

Re: ion addition by selective water replacement; addions.str
zee #23857 04/02/10 06:29 AM
Joined: Jan 2010
Posts: 15
Forum Member
Offline
Forum Member
Joined: Jan 2010
Posts: 15
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 -
Re: ion addition by selective water replacement; addions.str
Kellen #23863 04/02/10 06:33 PM
Joined: Sep 2003
Posts: 8,503
rmv Offline OP
Forum Member
OP Offline
Forum Member
Joined: Sep 2003
Posts: 8,503
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

Re: ion addition by selective water replacement; addions.str
rmv #25419 09/19/10 07:01 AM
Joined: Jul 2008
Posts: 12
T
Forum Member
Offline
Forum Member
T
Joined: Jul 2008
Posts: 12
Thanks for sharing it.
It's awesome

Re: ion addition by selective water replacement; addions.str
Todd #27241 04/18/11 06:05 PM
Joined: Mar 2011
Posts: 42
I
Forum Member
Offline
Forum Member
I
Joined: Mar 2011
Posts: 42
Hi,

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

thanks.

Re: ion addition by selective water replacement; addions.str
indian_user #27246 04/18/11 07:01 PM
Joined: Sep 2003
Posts: 4,796
Likes: 2
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 4,796
Likes: 2
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
Page 1 of 2 1 2

Moderated by  lennart, rmv 

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

PHP: 5.6.33-0+deb8u1 Page Time: 0.008s Queries: 35 (0.002s) Memory: 0.9907 MB (Peak: 1.1244 MB) Data Comp: Off Server Time: 2020-10-24 03:22:30 UTC
Valid HTML 5 and Valid CSS