CHARMM Development Project Forums CHARMM Community Script Archive applying crystal symmetry

The following is an example of applying CRYStal symmetry to create a unit cell from an asymmetric unit. The example is for a Cambridge Crystal Database structure of maltose, with P21 symmetry. The script is run twice, the first time to identify the image operation for the desired transform, the second time to apply the operation to the sugar and it's companion water molecule via COOR OPERate.

* make unit cell via coor oper; must run twice
* first time to identify 0,0,0 transforms, second to apply them
*

set ccdb maltos11
stream rtfprm.str

! SUGAR TOPOLOGY
read sequ card
* sugars
*
2
AGLC BGLC
gener g setup warn
patch 14ab g 1 g 2 setup
autogen dihe angl
ic purge
! JUST ADD WATER
read sequ tip3 1
gener w noang nodihe
join g w renum

! 3 LETTER PDB LIMITATION
rename resn glc sele resn aglc end
rename resn glc sele resn bglc end
rename resn tip sele resn tip3 end
read coor pdb name chm-@CCDB.pdb
rename resn aglc sele resid 1 end
rename resn bglc sele resid 2 end
rename resn tip3 sele resid 3 end

! CRYST1 4.866 15.077 10.702 90.00 97.07 90.00 P 21
set a 4.866
set b 15.077
set c 10.702
set beta 97.07
crystal define monoclinic @A @B @C 90. @BETA 90.
crystal build cutoff 12.0 noper 1
(-X,Y+1/2,-Z)

! IMAGE OP IS BASED ON TRANSFORM NUMBER FROM ABOVE LIST, e.g. C009
! UNCOMMENT STOP COMMAND FOR 1ST PASS, REPLACE COMMENT CHAR FOR 2ND
stop

! SHAKE X-H
shake bonh param
coor copy comp
cons fix sele .not. type H* end
coor shake
cons fix sele none end

ic fill
ic print

! REPEAT WITH LOWER SYMMETRY; GENERATE FULL UNIT CELL PSF, SAVE
gener t dupl g setup
join g t renum

open unit 3 write card name cell.psf
write psf unit 3 card
* @CCDB psf for unit cell
*
! DUPLICATE COORDS
coor dupl sele resid 1 end sele resid 4 end
coor dupl sele resid 2 end sele resid 5 end
coor dupl sele resid 3 end sele resid 6 end
! EDIT COOR OPER COMMANDS BASED ON OUTPUT FROM 1ST PASS; 0,0,0 TRANSFORMS
! CLOSEST AND SIMPLEST
coor oper C060 sele resid 4:6 end
! CENTER AT ORIGIN
coor orient norot

open unit 2 write card name cell0.crd
write coor card unit 2
* unit cell @CCDB initial construct
*
open unit 2 write card name cell0.pdb
write coor pdb unit 2
* unit cell @CCDB initial construct
*

! REBUILD WITH UNIT CELL; MUST BE SAME AS ABOVE ...
crystal free
crystal define monoclinic @A @B @C 90. @BETA 90.
crystal build cutoff 16.0 noper 0
! (-X,Y+1/2,-Z)

pressure init

coor copy comp
! BIAS TOWARD XTAL WITH RESTRAINTS
define asymhvy sele segid G .and. .not. type H* end
cons harm abso mass force 1.0 comp sele asymhvy end
mini abnr nstep 200 inbfrq 20 nprint 10 imgfrq 20 cutim 16. cutnb 16. -
vfswitch ctofnb 12.0 ctonnb 8. ewald pme order 6 spline kappa 0.34 -
fftx 6 ffty 16 fftz 10
cons harm abso mass force 0.5 comp sele asymhvy end
mini abnr nstep 100
cons harm clear
mini abnr nstep 100

pressure inst volume ?VOLU temp 295.

! MEASURE AND CALC DENSITY
coor stat
scalar mass stat
calc rho ( 1.66056 * ?STOT ) / ?VOLU

! COMPUTE RMSD FROM INITIAL COORDS
open read card unit 3 name cell0.crd
read coord card unit 3 comp
close unit 3
coor orient rms sele .not. type H* end comp

! OUTPUT RESULTING COORDS WITH SOME STATE INFO
open unit 1 write card name cellm.crd
write coor card unit 1
* unit cell @CCDB min e = ?ENER
* lattice min; rmsd = ?RMS ; p = ?PRSI ; density = @RHO
*
open unit 1 write card name cellm.pdb
write coor pdb unit 1
* unit cell @CCDB min e = ?ENER
* lattice min; rmsd = ?RMS ; p = ?PRSI ; density = @RHO
*

ic fill
open unit 3 write card name cellm.ic
write ic unit 3 card
* ic table for min unit cell @CCDB
*

stop

Dr. Venable,

I'm sorry if this is another newbie question. I'm still confuse about the way CHARMM handles periodic boundary condition. Lets say if I want to simulate only water molecule, is the CRYSTAL facility the only way to build the system with the images? If I don't use CRYSTAL facility, is there another way to do this?

There are two other ways beside CRYSTAL, the original IMAGE facility (which is used by CRYSTAL), and and a minimum image method known as BOUND; see images.doc

The CRYSTAL facility is required for the use of either P-M Ewald (EWALD PME) or NPT simulation methods.

Unless you have specific questions about the above script, posts in the topic-related forums are preferred.

Sorry about that. If I have more question definitely will ask in the proper forum. Thank you for the clarification about the CRYSTAL facility.

Hello Rick
I know the post was a long time ago. But I wonder why did you define crystal twice, first time with p21 space group and second time with P1 space group?

P21 is used both times. The second usage is because the PSF was changed, which invalidates the image data arrays.

Note that this script serves two purposes, and must be run twice. The COOR OPER command in the lower section must be edited based on the output listing from the first pass, which helps identify the transformation(s) needed to pack the unit cell from the asymmetric unit.

If I understand you correctly, once a unit cell is built, noper option is set to 0, since symmetry within the cell is already defined and noper=0 simply repeats the unit cell.

Yes, for the unit cell only the lattice constants are needed, so it is translation only (P1).