3.3. Density modification procedures

3.3.1. Map skeletonization

A get_skel.cmds command file:

MAP_DENS is the map that will be skeletonised, MAP_SKEL is a
connectivity map helping to create the skeleton atoms.  SKEL_CUT is
the level of density at which peaks and saddle points will be searched
according to S. Swanson's idea though the approach is somewhat
different and much faster.

> set vari MAP_DENS = 1
> set vari MAP_SKEL = 4
> set vari SKEL_CUT = 0.9

First remove old skeleton atoms and all symmetry equivalent segments.

> delete atom sele skel .or. segm name #* end
> key old sele all end

Find the extreme points in density.  A point is an EXTREME point if
all 6 neighboring points have a lower density value than the central
point.

> make map MAP_SKEL from MAP_DENS init 0 integer
> make point init from map MAP_DENS extreme SKEL_CUT 100.
> make map MAP_SKEL from MAP_DENS local SKEL_CUT 100.

Find clusters (clouds) of points around each EXTREME point and SADDLE
points between a two XTREME points.

> make map MAP_SKEL cloud

> make point from map MAP_DENS saddle MAP_SKEL

Create atoms and organize them in segments.  A segment is an entity
of skeleton atoms connected through network of bonds (extreme saddle
point connections).  All EXTREMEs that have no attached
skeleton atoms are DELETED.

> make atom from point reindex
> key skel sele .not old end
> dele atom sele .not by bond all .and skel end
> make segm from atom sele skel end

All bonds across the cells are deleted.  The Swanson's procedure is not
capable to link extreme points across the unit cell borders, but here in
MAIN segments are already organized according to their size, so we remove
here the cross-cell links to generate a nicer image.

> dele bond sele bond range 10. 1000. end
> dele atom sele .not old .a segm range 1 15 end
> set col sele skel .a segm range 10 100 end col 200
> set col sele skel .a segm range 100 10000 end col 240
> image set sele skel end set bond

> rena sele skel end atom CX
> set temp sele skel end 50.
> set weigh sele skel end 1.

> return

3.3.1. Solvent flattening - FFT procedure

Read an MIR map, symmetry operations and a reflection file.

> read file ../model/mir_clear_4.0_heavy.xmap map xplor
> read file >symm/c2.symm symm
> read file catl.fobs refl init limit -99 0 0 reso 3.5 100.0

Set variables MAP_DENS (origin of density) and MAP_SOLV (solvent
flattened map).

> set vari MAP_DENS = 1
> set vari MAP_SOLV = 2
> set vari SOLV_CUT = 0.4
> set vari SOLV_RAD = 10.0
> set vari HOLE_SIZE = 6

> make map MAP_SOLV from MAP_DENS init -9999  cell real

The second command line invokes the usual real space procedure,
which is maybe even 100 times slower that the one described below.

> make map MAP_SOLV from MAP_DENS init -9999  cell real
> make map MAP_SOLV from MAP_DENS cut 0.0 radi SOLV_RAD weight

Initiate a unit cell with the density set to EMPTY and a density of an
atom that is placed at coordinate system origin and which density is
distance dependent (1. - d / SOLV_RAD).  Set the empty points to zero,
Fourier transform the map and save its Fourier coeficients in
the FWORKSET array.

> make map MAP_SOLV from MAP_DENS init -9999  cell real
> make map MAP_SOLV atom dist SOLV_RAD
> make map MAP_SOLV set -100000 0 0.0
> fourier map MAP_SOLV
> reflect set ampl phase fwork = fcalc * 1.0

Copy the MIR density into the new MAP_SOLV, set all negative density
values to zero and Fourier transform it.  R-value calculation is
performed in order to follow convergence of a further procedure.

> make map MAP_SOLV from 2 init 9999  cell real
> make map MAP_SOLV from MAP_DENS copy
> make map MAP_SOLV set -99999 0 0.0
> fourier map MAP_SOLV
> reflect shells 10 r-values 

Multiply FWORKSET (sphere model) and FCALCULATE (truncated MIR
density) Fourier coeficients, convert the MAP_SOLV map to a convex
map, fill it with the multiplied Fourier coefficients and Fourier BACK
transform it.

> reflect set ampl phase fwork = fwork * 1.0  *  fcalc * 1.0
> make map MAP_SOLV conv complex
> make map MAP_SOLV zero
> refl fill-map MAP_SOLV defined
> four map MAP_SOLV back 
> make map MAP_SOLV rescale 

From here on the real space and FFT driven procedures continue
together.  Perform the score map analyses (map occupancy in 256
steps), convert the selected parts to a MASK region and copy density
into it.  Set the density of the empty region to, Fourier transform it
and caluclate R-values.

> anal ini
> anal map MAP_SOLV
> make map MAP_SOLV anal range SOLV_CUT

Here we can in addition to range play also with the enevelope cavity size:

> make map MAP_SOLV + 1 from MAP_SOLV init 9999 copy
> make map MAP_SOLV from MAP_SOLV + 1 fill HOLE_SIZE


> make map MAP_SOLV from MAP_DENS copy
> make map MAP_SOLV set -99999 -9000 0.0
> fourier map MAP_SOLV
> reflect shells 10 r-values 

Save the FCALCULATE values in the FWORKSET array and calculate
an Fobs map 

> reflect set ampl phase fwork = fcalc * 1.0
> reflect set ampl fwork = fobs * 1.0

or 2 Fobs - Fcalc map.

> reflect set ampl scale fcalc fwork = fobs * 2.0 - fcalc * 1.0


> make map MAP_SOLV set zero
> refl fill-map MAP_SOLV defined 
> fourier back map MAP_SOLV
> make map MAP_SOLV rescale

> return

Here you can either save or display the solvent flattened
map or calculated the next cycle by redefining the variables:

> set vari MAP_DENS = 2
> set vari MAP_SOLV = 3


The menu driven solvent flattenig procedure command file
solv_flat.cmds and solv_flat_fast.com (essentially described above):

! Set variables MAP_DENS (origin of density) and MAP_SOLV (solvent
! falttened map).

> set vari MAP_DENS = 2
> set vari MAP_SOLV = 3
> set vari SOLV_CUT = 0.4
> set vari SOLV_RAD = 8.0
> <solv_flat_fast
> return

Redefine variables to perform the next cycle.

> set vari MAP_DENS = 2
> set vari MAP_SOLV = 3
> <solv_flat_fast
> return

3.3. Electron density averaging

See the main/examples/aver_n_mol directory and the MAIN manual
where the procedure is described in details.  

3.4. Atom density correlation

You can essentially take the solv_flat_fast procedure, only 
the command line

> make map MAP_SOLV atom dist SOLV_RAD

should be replaced with

> read coor
> *
> CORI

> key origin sele atom numb natom end
> set weigh sele origin end 1.0 ! or 5.0 if you wish
> set temp sele origin end 40. ! or some other values 
> make map MAP_SOLV atom dens sele atom numb natom end \
>                   atom symmetry 1

which creates an atom at the origin of coordinate system.

If you want to calculate your model or density translation function
you can replace these lines by reading your molecular model
or instead of filling the MAP_SOLV map with atomic density
offer electron density map or its part.