atom

Atomistic Topology Operations in MATLAB (atom), is a MATLAB library for manipulation of (periodic) molecular systems

Atomistic Topology Operations in MATLAB, functions for manipulation of molecular dynamics or monte carlo simulation systems.

% Note that version 2.11 comes with an extensive and searchable html-documentation for all the >100 functions, which can be used interactively from MATLAB's own browser.

% The purpose of the atom library is to automate and enable efficient construction/manipulation and analysis of complex and multicomponent molecular systems, and generate topological information with bonds and angles etc, over the PBC!

It is especially useful for building inorganic/geochemical systems, since bond distances can be compared to the ideal semi-empirical bond distances computed with the Bond Valence Sum Method, or just simply just compared to Shannon's revised radii. Or one could plot a theoretical X-ray profile with the xrd_atom() function.

% For lists of all available functions by category, see inside these files:

List_all_functions.m

List_build_functions.m

List_export_functions.m

List_general_functions.m

List_import_functions.m

List_forcefield_functions.m

% The atom scripts can read and write basic .pdb | .xyz | .gro | .mol2 structure files as well as write basic .itp and .psf topology files. They can also manipulate/transform the structures in various ways making use of the MATLAB struct variable and indexing. They can be used to build and plot multicomponent systems, by adding molecules, ions and SPC/TIP3P/TIP4P water molecules or other solvents (ie solvating an existing molecule/slab) into a simulation box, and remove molecular overlap. For plotting one can call vmd(atom,Box_dim) if the VMD software is also installed and the PATH2VMD() function is properly set. Else the very quick-and-dirty plot_atom(arguments) or the slower show_atom(arguments) can be used. Most functions takes PBC into account, which allows for generation of topologies of molecules with bonds, angles, dihedrals across the PBC. There is also some support for triclinic support using the tilt vectors xy, xz, yz, as well as for generating powder X-ray diffractograms using the function xrd_atom(). The Bond Valence Sum method is also implemented in the analyze_atom() function.

Michael Holmboe 

[e-mail] 

Chemistry department

Umeå University, Sweden

% Where to start? Perhaps the html-documentation with some basic examples? 

% 

% Some typical commands...

%

% To read a structure file into matlab (check the variable explorer) 

atom=import_atom('filename') % filename could be a .pdb | .xyz | .gro file 

% or... 

atom=import_atom_pdb('filenamepdb')

atom=import_atom_gro('filenamegro') 

atom=import_atom_xyz('filenamexyz')

% Note that you get a lot more info then just the atom struct variable, like the box dimension variable Box_dim 

% To write a atom struct to a new topology or structure file

write_atom_lmp(atom,Box_dim,filename,1.2,1.2,'clayff') % supports bonds, angles, simple dihedrals

write_atom_psf(atom,Box_dim,filename,1.2,1.2,'clayff') % note only bonds and angles

write_atom_itp(atom,Box_dim,filename,1.2,1.2,'clayff','spce') % Gromacs topology file, note only bonds and angles

write_atom_pdb(atom,Box_dim,filename)

write_atom_cif(atom,Box_dim,filename)

write_atom_gro(atom,Box_dim,filename) 

write_atom_xyz(atom,Box_dim,filename)

% Adding water to a box 

% - This function SOLvates a certain region defined by limits with a water or a custom solvent

% structure with density. r (and r-0.5 for H) is the closest distance of solvent atoms

% to the (optional) solute atoms

SOL_atom = solvate_atom(limits,density,r,maxsol) % limits can be [10] | [10 20 30] | [10 20 30 40 50 60]

SOL_atom = solvate_atom(limits,density,r,maxsol,solute_atom,'tip4p') % spc | tip3p | tip4p | tip5p

% One can filter the atom struct with respect to molid, resname, atomtype, index, coordinates and so on. This allows manipulation of an atom struct on the atomic, molecule and molecular type level. This also allows us to use  'dynamic indexes' of groups of atom.{molid/resname/type/index/} when analyzing a trajectory for instance. Some basic examples: 

index=ismember([atom.type],[{'Al' 'Alt' 'Mgo'}]) % gives a binary (1/0) logical array 

index=strcmp([atom.type],'Al') % try also strncmp or strncmpi? 

index=find(strncmpi([atom.type],'al',2) % Will find the indexes of 'Al' 'Alt? 

new_atom=atom(index) % This creates a new_atom struct with the filtered/selected atomtypes 

positive_z_atom=atom([atom.z]>0) % finds all atoms with a positve z-coordinate

first100_atom=atom([atom.index]<101) % finds the first 100 atoms in the atom struct 

first100_v2_atom=atom(1:100) % also finds the first 100 atoms in the atom struct

% Merging two different atom structs 

% - This function returns the second atom set with non-overlapping atoms

new_atom = merge_atom(atom1,Box1,atom2,Box2,type,Atom_label,r)

% Calculating bonds or the distance matrix/es 

atom = bond_angle_dihedral_atom(atom,Box_dim,1.2,2.2) % first cut-off for bonds to H's, second cut-off for everything else. Only atoms having the same MolID will be considered as bonded.

dist_matrix = dist_matrix_atom(atom,Box_dim) % another cell-list-version also exist for orthogonal systems.