CHARMM Development Project Forums CHARMM Community Script Archive ion addition by selective water replacement

Dear Users,

At first, thank You this really useful script!

I have modifed this 'addions.str' script, because I would like to replace the water molecules with NO gas molecules. Now it looks like this:

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

! PICK A WATER MOLECULE AND REPLACE IT WITH A CHLORIDE ION
set k 1
label ilp
define xcld sele ( @SOL ) .or. segid NO end ! EXCLUDED SEGMENTS
define prox sele xcld .around. @MND end ! NEARBY EXCLUDED ATOMS
define list sele atom SOLV * 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 NO @K NO end
delete atom sele .byres. targ end sort ! REMOVE THE WATER MOLECULE
incr k by 1
if k .le. @ICL goto ilp


return

When I try with this one, a warning level occurs in the 'loop*.log' file:

CHARMM> coor set xdir ?XAVE ydir ?YAVE zdir ?ZAVE sele atom NO @K NO end
Parameter: K -> "1"
RDCMND substituted energy or value "?XAVE" to "33.118"
RDCMND substituted energy or value "?YAVE" to "6.324"
RDCMND substituted energy or value "?ZAVE" to "17.03"
SELRPN> 0 atoms have been selected out of 32947

***** LEVEL 0 WARNING FROM *****
***** ZERO ATOMS SELECTED
******************************************
BOMLEV ( -1) IS NOT REACHED. WRNLEV IS 5

SELECTED COORDINATES SET IN THE MAIN SET.

My NO residue is a three-site gas model:
RESI NO 0.000
GROUP
ATOM OM LP 0.595
ATOM N NO -0.250
ATOM O ON -0.345
BOND N O
ACCEPTOR O
ACCEPTOR N
PATCHING FIRS NONE LAST NONE


Can anyone explains me why CHARMM can not do the selection?

Thank You in advance,
Aniko

There is no atom named NO in your NO residue. The atom names are OM, N and O.

The COOR SET command is intended for a single atom, not a 3-atom molecule, so further changes will be needed. Both the initial COOR SET, which put ions at the origin, and the one in the loop need to be changed, probably to COOR TRANSlate.

For placing multi-atom residues, I would probably read the coordinates for each molecule from a file, and then translate it to the desired location. There are some nuances in that depending on the type of file used, and the initial coordinates of the molecule(s) being read.

Thank You Lennart!

I have replaced with the atom name of N. Now the script works. But because I replaced only one atom, the OH2, just for this N atom has coordinates, the O and OM have not, of course. I would like to assign coordinates for O and OM too. If I replace the one of the H atoms of the water molecules by O, there are also no coordinates for OM. Anyway how can a I do this? How can I define a second atom selection?

I have also tried to write out pdb file from this running and regenerate the structure from a pdb which contains just the N atoms. In that case, CHARMM added 9999.0009999.0009999.000 coordinates for the O and OM.

Do You have any suggestion have can I solve this problem?

Thank You in advance,
Aniko

Dear Rick,

So if I understand well, the best solution is if I generate a pdb file which contains 1 NO molecule with coordinates and than I should translate them in every loop cycle. In order to do this I have modified the addions.str in the following way:

read sequ no @ICL
generate no noang nodihe first none last none
coor set xdir -0.031 ydir 4.923 zdir 17.073 sele atom NO * OM end
coor set xdir -0.183 ydir 5.621 zdir 16.975 sele atom NO * N end
coor set xdir 0.064 ydir 4.485 zdir 17.134 sele atom NO * O end

! PICK A WATER MOLECULE AND REPLACE IT WITH A CHLORIDE ION
set k 1
label ilp
define xcld sele ( @SOL ) .or. segid NO end ! EXCLUDED SEGMENTS
define prox sele xcld .around. @MND end ! NEARBY EXCLUDED ATOMS
define list sele atom SOLV * 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 NO @K N end
coor trans xdir ?XAVE ydir ?YAVE zdir ?ZAVE sele atom NO @K OM end
coor trans xdir ?XAVE ydir ?YAVE zdir ?ZAVE sele atom NO @K N end
coor trans xdir ?XAVE ydir ?YAVE zdir ?ZAVE sele atom NO @K O end
delete atom sele .byres. targ end sort ! REMOVE THE WATER MOLECULE
incr k by 1
if k .le. @ICL goto ilp


return


Now this script works, replaces the 188 NO molecules with real coordinates of each of the atoms, but several of this newly built molecules are out of the water box. Do You have any suggestion why happens this?

Thank You in advance,
Aniko

Your template NO molecules are not placed at the origin of the coordinate system, so your translate commands don't put them where the water molecule was.
Try adding this line before !PICK A WATER ...
coor trans xdir -0.064 ydir -4.485 zdir -17.134 sele segi NO end

Dear Lennart (and Rick),

Now it works well!!! Thank You very much!!!

Best wishes,
Aniko