Topic Options
#3455 - 10/01/04 11:25 AM Modified solvent-box.str
davhak Offline
Forum Member

Registered: 07/05/04
Posts: 162
Loc: 40.1N, 44.3E
Dear All,

Below is the modified solvent-box.str script, which automatically creates new segment names (starting with "WT") to avoid the residue excess per segment (over 9999). For sake of readability the old code is left commented and changes are placed right after the commented code lines.

Code:


*NAME: solvent-box.str
*PURPOSE: make a TIP3 box with user-defined edges, and make room for solute
*AUTHOR: Lennart Nilsson, Karolinska Institutet (October 7, 2003)
* THERE IS NO ERROR CHECKING IN THIS FILE
* Dynamic segment generation added for those CHARMM versions, which do not
* support more than 9999 residues per segment . Modification
* by Davit Hakobyan (October 1, 2004).
*
!ASSUMPTIONS:
! The psf and coordinates of the solute are present
! A selection "SOLUTE" defining the solute has been made
! (define solute sele ... end)
! There is no segment present with segid WAT, WAT0, WAT1, WAT2 or WAT3
! XSIZ, YSIZ and ZSIZ have been s set to the desired box dimensions
! Uses $CHM_STREAM/read-watbox.str to get coordinates for primary water box
! Unix environment variable CHM_STREAM points to directory containing these
! scripts
!RESULT:
! At return there will be a new segment named WAT containing TIP3 water
! molecules ! that have their oxygen atom in the box, and not overlapping
! with the solute.
! The system is centered at the origin, with the main axis of inertia of the
! SOLUTE! along the x-axis. For the deletion of overlapping waters only
! heavy solute atoms and water oxygens are considered; overlap distance
! criteria are defined below
! No other manipulations (minimizations etc) are performed on the system;
! bomblevel is set to -2 to avoid stopping at empty selections or
! empty segments
! 2004.10.01 modification
! For those segments, which threaten to contain more than 9999 residues per segment
! new segment name is created. This results in some changes in water segments
! names. For current modification water segments start with name WT followed by
! a number starting from 2. Due to intermediate deletion and generation of atoms
! and segments the final structure may not contain water segment starting from
! strictly WT2 and may contain gap in segment naming. Modification could be mostly
! found if search for "Auto_naming" comment line.

!USAGE EXAMPLE:
!read rtf card name /applic/charmm/toppar/top_all22_prot.inp
!read para card name /applic/charmm/toppar/par_all22_prot.inp
!read psf card name my_structure.psf
!read coor card name my_structure.crd
!set XSIZ 27.0
!set YSIZ 30.0
!set ZSIZ 35.0
!define SOLUTE sele segid prot .or. segid dna1 .or. segid dna2 end
!stream $CHM_STREAM/solvent-box.str
!! here you may want to add some boundary potential
!cons harm force 50.0 sele SOLUTE end
!mini sd nstep 50
!mini abnr nstep 50
!. . .
!LOCAL VARIABLES:
! XTR,YTR,ZTR,XNUM,YNUM,ZNUM,BOX,XB2,YB2,ZB2
! NW0,NROD,NSLAB,OVERLAP,OVERAP2,I
!Defines selection SOLUTE1
! distance to use when deleting waters that overlap with the solute
! possibly different for second overlay of slightly shifted/rotated slab
set OVERLAP 2.2
set OVERLAP2 2.2
! solvation will be done by sliding a slab of TIP3 along the x-axis
! the water box we start with is a cube with 18.856A side containing 216 TIP3
set BOX 18.856
set NW0 216

! Auto_naming
set NLIMIT 9999 ! Maximum number of residues per segment
set swi 2 ! Starting Water Index constant
set wiCur @swi ! Variable which keeps the current index of water segment
! along this script.

! for the water-deletions we will use only the heavy atoms of the solute
define SOLUTE1 sele SOLUTE .and. .not. hydrogen end
! set bomblevel to -2 to avoid stopping on empty selections or empty segments
bomb -2

! begin calculations
! place SOLUTE along xaxis centered at (XSIZ/2,YSIZ/2,ZSIZ/2) in first octant so
! we can easily move our box around to cover the solute
calc XB2 = @XSIZ / 2.0
calc YB2 = @YSIZ / 2.0
calc ZB2 = @ZSIZ / 2.0
coor orient sele SOLUTE end
coor trans xdir @XB2 ydir @YB2 zdir @ZB2

calc XNUM = INT(@XSIZ / @BOX + 1)
calc YNUM = INT(@YSIZ / @BOX + 1)
calc ZNUM = INT(@ZSIZ / @BOX + 1)

! number of water molecules in a rod and in a slab
calc NROD = @NW0 * @YNUM
calc NSLAB = @ZNUM * @NROD

!generate one rod of boxes along Y-axis
read sequence tip3 216
gene wat0 noangle nodihe
stream $CHM_STREAM/read-watbox.str
! put it in the first octant
coor stat sele segid wat0 .and. .not. hydrogen end
calc XTR = - ?XMIN
calc YTR = - ?YMIN
calc ZTR = - ?ZMIN
coor trans xdir @XTR ydir @YTR zdir @ZTR sele segid WAT0 end

! Auto_naming
read sequence tip3 216
!gene wat noangle nodihe
gene wt@wiCur noangle nodihe
!coor dupl sele segid wat0 end sele segid wat end
coor dupl sele segid wat0 end sele segid wt@wiCur end

if YNUM .eq. 1 goto do_z
! loop over the remaining 2 to NUM boxes in the rod
set I 2

! Auto_naming
set wiSum @NW0

label yloop
read sequence tip3 216
gene wat1 noangle nodihe
coor dupl sele segid wat0 end sele segid wat1 end
coor transl ydir @BOX sele segid wat1 end
coor transl ydir @BOX sele segid wat0 end
!join wat wat1 renumber

! Auto_naming
if @NROD .le. @NLIMIT THEN GOTO JUNC_PNT1
set wiLast @wiCur
calc wiSum = @wiSum + @NW0
calc wiCur = INT(@wiSum / @NLIMIT + @wiCur)
if @wiLast .eq. @wiCur then goto JUNC_PNT1
rename segid wt@wiCur sele segid wat1 end
label JUNC_PNT1
join wt@wiCur wat1 renumber

incr I by 1
if I le @YNUM goto yloop

! done with the rod, make a slab perpendicular to X-axis
label DO_Z
! wat0 is no longer needed
delete atom sele segid wat0 end
if ZNUM .eq. 1 goto do_x

!rename segid wat0 sele segid wat end

!read sequence tip3 @NROD
!gene wat noangle nodihe
!coor dupl sele segid wat0 end sele segid wat end

! Auto_naming
set wiLast @wiCur
set wiSum @NROD
set wiTmp @swi
calc NumTimes = INT(@NLIMIT / @NW0)
calc AtmBlock = @NumTimes * @NW0

label GEN_PNT1
if wiSum .le. 0 then goto PRE_ZLOOP
if wiSum .le. @AtmBlock then goto GEN_PNT2
set AtmNum @AtmBlock
goto GEN_PNT3
label GEN_PNT2
set AtmNum @wiSum
label GEN_PNT3
increment wiCur by 1
read sequence tip3 @AtmNum
gene wt@wiCur noangle nodihe
coor dupl sele segid wt@wiTmp end sele segid wt@wiCur end
decrement wiSum by @AtmNum
increment wiTmp by 1
goto GEN_PNT1

label PRE_ZLOOP

set I 2
label ZLOOP
!read sequence tip3 @NROD
!gene wat1 noangle nodihe
!coor dupl sele segid wat0 end sele segid wat1 end
!coor transl zdir @BOX sele segid wat1 end
!coor transl zdir @BOX sele segid wat0 end
!join wat wat1 renumber

! Auto_naming
set wiSum @NROD
set wiTmp @swi
label GEN_PNT10
if wiSum .le. 0 then goto JUNC_PNT2
if wiSum .le. @AtmBlock then goto GEN_PNT20
set AtmNum @AtmBlock
goto GEN_PNT30
label GEN_PNT20
set AtmNum @wiSum
label GEN_PNT30
increment wiCur by 1
read sequence tip3 @AtmNum
gene wt@wiCur noangle nodihe
coor dupl sele segid wt@wiTmp end sele segid wt@wiCur end
coor transl zdir @BOX sele segid wt@wiCur end
coor transl zdir @BOX sele segid wt@wiTmp end
decrement wiSum by @AtmNum
increment wiTmp by 1
goto GEN_PNT10

label JUNC_PNT2
incr I by 1
if I le @ZNUM goto zloop

! done with the slab, now let it slide along the x-axis, after centering everything at orig
label DO_X
!delete atom sele segid wat0 end

! Auto_naming
set wiTmp @swi
label DEL_PNT1
delete atom sele segid wt@wiTmp end
increment wiTmp by 1
if wiTmp .le. @wiLast then goto DEL_PNT1
calc wiFirst = @wiLast + 1

!rename segid wat0 sele segid wat end

coor transl xdir -@XB2 ydir -@YB2 zdir -@ZB2

! Auto_naming
! This auto naming snippet performs packing of those water segments
! which might fit several of others. Also renames the names of the
! final compact segment to start from the starting index (defined
! by "swi" variable
calc NumTimes = INT(@NLIMIT / @NROD)
set wiTmp @swi
label SHIFT_PNT1
set wFit 2
calc wiNext = @wiFirst + 1
if wiFirst .ge. @wiCur then goto SHIFT_PNT4
label SHIFT_PNT2
if wFit .gt. @NumTimes then goto SHIFT_PNT3
join wt@wiFirst wt@wiNext
increment wiNext by 1
increment wFit by 1
goto SHIFT_PNT2
label SHIFT_PNT3
rename segid wt@wiTmp sele segid wt@wiFirst end
increment wiTmp by 1
set wiFirst @wiNext
goto SHIFT_PNT1
label SHIFT_PNT4
rename segid wt@wiTmp sele segid wt@wiFirst end
set wiCur @wiTmp
set wiLast @wiCur

!read sequence tip3 @NSLAB
!gene wat2 noangle nodihe
!coor dupl sele segid wat0 end sele segid wat2 end
! delete water molecules further away than XB2, YB2, ZB2 A from origin
!delete atom sele .byres. ( segid wat2 .and. type OH2 .and. -
! ( property abs x .gt. @XB2 .or. property abs y .gt. @YB2 -
! .or. property abs z .gt. @ZB2) ) end
! and also those overlapping with the SOLUTE
!delete atom sele .byres. ( segid wat2 .and. type OH2 .and. -
! ( SOLUTE1 .around. @OVERLAP ) ) end

! Auto_naming
set wiSum @NSLAB
set wiTmp @swi
calc AtmBlock = @NumTimes * @NROD
label GEN_PNT100
if wiSum .le. 0 then goto DEL_PNT2
if wiSum .le. @AtmBlock then goto GEN_PNT200
set AtmNum @AtmBlock
goto GEN_PNT300
label GEN_PNT200
set AtmNum @wiSum
label GEN_PNT300
increment wiCur by 1
read sequence tip3 @AtmNum
gene wt@wiCur noangle nodihe
coor dupl sele segid wt@wiTmp end sele segid wt@wiCur end
! delete water molecules further away than XB2, YB2, ZB2 A from origin
delete atom sele .byres. ( segid wt@wiCur .and. type OH2 .and. -
( property abs x .gt. @XB2 .or. property abs y .gt. @YB2 -
.or. property abs z .gt. @ZB2) ) end
! and also those overlapping with the SOLUTE
delete atom sele .byres. ( segid wt@wiCur .and. type OH2 .and. -
( SOLUTE1 .around. @OVERLAP ) ) end
decrement wiSum by @AtmNum
increment wiTmp by 1
goto GEN_PNT100

label DEL_PNT2
! do the same slab again, but translated -1A and rotated 30 degrees around the x-axis
!read sequence tip3 @NSLAB
!gene wat3 noangle nodihe
!coor dupl sele segid wat0 end sele segid wat3 end
!coor rota xdir 1.0 phi 30.0 sele segid wat3 end
!coor transl xdir -1.0 sele segid wat3 end
! delete water molecules further away than XB2, YB2, ZB2 A from origin
!delete atom sele .byres. ( segid wat3 .and. type OH2 .and. -
! ( property abs x .gt. @XB2 .or. property abs y .gt. @YB2 -
! .or. property abs z .gt. @ZB2) ) end
! and also those overlapping with the SOLUTE or existing WAT2 segment
!delete atom sele .byres. ( segid wat3 .and. type OH2 .and. -
!( (SOLUTE1 .or. (segid wat2 .and. type OH2)) .around. @OVERLAP2 ) ) end
!join wat2 wat3 renumber

! Auto_naming
set wiSum @NSLAB
set wiTmp @swi
label GEN_PNT400
if wiSum .le. 0 then goto PRE_XLOOP
if wiSum .le. @AtmBlock then goto GEN_PNT500
set AtmNum @AtmBlock
goto GEN_PNT600
label GEN_PNT500
set AtmNum @wiSum
label GEN_PNT600
increment wiCur by 1
read sequence tip3 @AtmNum
gene wt@wiCur noangle nodihe
coor dupl sele segid wt@wiTmp end sele segid wt@wiCur end
coor rota xdir 1.0 phi 30.0 sele segid wt@wiCur end
coor transl xdir -1.0 sele segid wt@wiCur end
! delete water molecules further away than XB2, YB2, ZB2 A from origin
delete atom sele .byres. ( segid wt@wiCur .and. type OH2 .and. -
( property abs x .gt. @XB2 .or. property abs y .gt. @YB2 -
.or. property abs z .gt. @ZB2) ) end
! and also those overlapping with the SOLUTE or existing WAT2 segment
calc wiFirst = @wiLast + 1
label DEL_PNT4
delete atom sele .byres. ( segid wt@wiCur .and. type OH2 .and. -
( (SOLUTE1 .or. (segid wt@wiFirst .and. type OH2)) .around. @OVERLAP2 ) ) end
increment wiFirst by 1
if wiFirst .lt. @wiCur then goto DEL_PNT4

decrement wiSum by @AtmNum
increment wiTmp by 1
goto GEN_PNT400

label PRE_XLOOP

if XNUM .eq. 1 goto done

set I 2
label XLOOP
! read sequence tip3 @NSLAB
! gene wat3 noangle nodihe
! coor dupl sele segid wat0 end sele segid wat3 end
! coor transl xdir @BOX sele segid wat3 end
! coor transl xdir @BOX sele segid wat0 end
! delete water molecules further away than XB2, YB2, ZB2 A from origin
! delete atom sele .byres. ( segid wat3 .and. type OH2 .and. -
! ( property abs x .gt. @XB2 .or. property abs y .gt. @YB2 -
! .or. property abs z .gt. @ZB2) ) end
! and also those overlapping with the SOLUTE
! delete atom sele .byres. ( segid wat3 .and. type OH2 .and. -
! ( SOLUTE1 .around. @OVERLAP ) ) end
! join wat2 wat3 renumber

! Auto_naming
set wiSum @NSLAB
set wiTmp @swi
label GEN_PNT700
if wiSum .le. 0 then goto GEN_PNT1000
if wiSum .le. @AtmBlock then goto GEN_PNT800
set AtmNum @AtmBlock
goto GEN_PNT900
label GEN_PNT800
set AtmNum @wiSum
label GEN_PNT900
increment wiCur by 1
read sequence tip3 @AtmNum
gene wt@wiCur noangle nodihe
coor dupl sele segid wt@wiTmp end sele segid wt@wiCur end
coor transl xdir @BOX sele segid wt@wiCur end
coor transl xdir @BOX sele segid wt@wiTmp end
! delete water molecules further away than XB2, YB2, ZB2 A from origin
delete atom sele .byres. ( segid wt@wiCur .and. type OH2 .and. -
( property abs x .gt. @XB2 .or. property abs y .gt. @YB2 -
.or. property abs z .gt. @ZB2) ) end
! and also those overlapping with the SOLUTE
delete atom sele .byres. ( segid wt@wiCur .and. type OH2 .and. -
( SOLUTE1 .around. @OVERLAP ) ) end
! and also those overlapping with the SOLUTE
delete atom sele .byres. ( segid wt@wiCur .and. type OH2 .and. -
( SOLUTE1 .around. @OVERLAP ) ) end
decrement wiSum by @AtmNum
increment wiTmp by 1
goto GEN_PNT700

label GEN_PNT1000
! do the same slab again, but translated -1A and rotated 30 degrees around the x-axis
! read sequence tip3 @NSLAB
! gene wat3 noangle nodihe
! coor dupl sele segid wat0 end sele segid wat3 end
! coor rota xdir 1.0 phi 30.0 sele segid wat3 end
! coor transl xdir -1.0 sele segid wat3 end
! delete water molecules further away than XB2, YB2, ZB2 A from origin
! delete atom sele .byres. ( segid wat3 .and. type OH2 .and. -
! ( property abs x .gt. @XB2 .or. property abs y .gt. @YB2 -
! .or. property abs z .gt. @ZB2) ) end
! and also those overlapping with the SOLUTE or existing WAT2 segment
! delete atom sele .byres. ( segid wat3 .and. type OH2 .and. -
! ( (SOLUTE1 .or. (segid wat2 .and. type OH2) ).around. @OVERLAP2 ) ) end
! join wat2 wat3 renumber

! Auto_naming
set wiSum @NSLAB
set wiTmp @swi
label GEN_PNT1100
if wiSum .le. 0 then goto JUNC_PNT3
if wiSum .le. @AtmBlock then goto GEN_PNT1200
set AtmNum @AtmBlock
goto GEN_PNT1300
label GEN_PNT1200
set AtmNum @wiSum
label GEN_PNT1300
increment wiCur by 1
read sequence tip3 @AtmNum
gene wt@wiCur noangle nodihe
coor dupl sele segid wt@wiTmp end sele segid wt@wiCur end
coor rota xdir 1.0 phi 30.0 sele segid wt@wiCur end
coor transl xdir -1.0 sele segid wt@wiCur end
! delete water molecules further away than XB2, YB2, ZB2 A from origin
delete atom sele .byres. ( segid wt@wiCur .and. type OH2 .and. -
( property abs x .gt. @XB2 .or. property abs y .gt. @YB2 -
.or. property abs z .gt. @ZB2) ) end
! and also those overlapping with the SOLUTE or existing WAT2 segment
calc wiFirst = @wiLast + 1
label DEL_PNT5
delete atom sele .byres. ( segid wt@wiCur .and. type OH2 .and. -
( (SOLUTE1 .or. (segid wt@wiFirst .and. type OH2) ).around. @OVERLAP2 ) ) end
increment wiFirst by 1
if wiFirst .lt. @wiCur then goto DEL_PNT5

decrement wiSum by @AtmNum
increment wiTmp by 1
goto GEN_PNT1100

label JUNC_PNT3
incr I by 1
if I le @XNUM goto xloop

label DONE
!delete atom sele segid wat0 end
! Auto_naming
set wiTmp @swi
label DEL_PNT6
delete atom sele segid wt@wiTmp end
increment wiTmp by 1
if wiTmp .le. @wiLast then goto DEL_PNT6
!calc wiFirst = @wiLast + 1

!rename segid wat sele segid wat2 end
return



Top
#3456 - 12/08/04 01:09 AM Re: Modified solvent-box.str [Re: davhak]
chuan Offline
Forum Member

Registered: 05/17/04
Posts: 222
This script works perfectly.
_________________________
CHARMM 30b1 driven by 1/ Xeon (32 bits) 2/ Redhat 7.3 (32 bits) with a Quadrics-modified 2.4-18-5 kernel 3/ Chuan, with 95% of the mentorship coming from great scientists frequenting this forum. 4/ Gracious support from the forum.

Top

Moderator:  chmgr, John Legato, petrella