---

BaseMolecule.h

Go to the documentation of this file.

/***************************************************************************
 *cr                                                                       
 *cr            (C) Copyright 1995-2019 The Board of Trustees of the           
 *cr                        University of Illinois                       
 *cr                         All Rights Reserved                        
 *cr                                                                   
 ***************************************************************************/

/***************************************************************************
 * RCS INFORMATION:
 *
 *      $RCSfile: BaseMolecule.h,v $
 *      $Author: johns $        $Locker:  $             $State: Exp $
 *      $Revision: 1.151 $      $Date: 2020/10/28 15:09:56 $
 *
 ***************************************************************************/
#ifndef BASEMOLECULE_H
#define BASEMOLECULE_H

#ifndef NAMELIST_TEMPLATE_H
#include "NameList.h"
#endif
#ifndef RESIZEARRAY_TEMPLATE_H
#include "ResizeArray.h"
#endif
#include "Atom.h"
#include "Residue.h"
#include "Timestep.h"
#include "Fragment.h"
#include "Matrix4.h"
#include "intstack.h"

#ifdef VMDWITHCARBS
#include "SmallRing.h"
#include "SmallRingLinkages.h"
#include "inthash.h"
#endif

#if 0
#define VMDFASTRIBBONS
#endif

class VolumetricData;
class QMData;

class BaseMolecule {
public:
  //
  // public molecular structure data (for ease of access):
  //
  int nAtoms;                         
  int nResidues;                      
  int nWaters;                        
  int nSegments;                      
  int nFragments;                     
  int nProteinFragments;              
  int nNucleicFragments;              

  NameList<int> atomNames;            
  NameList<int> atomTypes;            
  NameList<int> resNames;             
  NameList<int> chainNames;           
  NameList<int> segNames;             
  NameList<int> altlocNames;          

  ResizeArray<Residue *> residueList; 
  ResizeArray<Fragment *> fragList;   

  ResizeArray<Fragment *> pfragList;  

  ResizeArray<int> pfragCyclic;       
#if defined(VMDFASTRIBBONS)
  ResizeArray<int *> pfragCPList;     

#endif


  ResizeArray<Fragment *> nfragList;  
  ResizeArray<int> nfragCyclic;       
#if defined(VMDFASTRIBBONS)
  ResizeArray<int *> nfragCPList;     

#endif

#ifdef VMDWITHCARBS
  ResizeArray<SmallRing *> smallringList; 

  SmallRingLinkages smallringLinkages; 
  int currentMaxRingSize;             
  int currentMaxPathLength;           
#endif

  // Molecule instance transformation matrices, not including self
  ResizeArray<Matrix4> instances;     

  // Extra molecular dynamics structure info not normally needed for
  // anything but structure building tasks
  ResizeArray<int> angles;            
  ResizeArray<int> dihedrals;         
  ResizeArray<int> impropers;         
  ResizeArray<int> cterms;            

  // Extra (optional) data for this molecule.
  // If the name is unassigned then the values can be taken to be zero.
  NameList<float *> extraflt;         
  NameList<int *> extraint;           
  NameList<unsigned char *> extraflg; 

  // more structure topology handling info. assign names to FF types.
  ResizeArray<int> angleTypes;        
  ResizeArray<int> dihedralTypes;     
  ResizeArray<int> improperTypes;     
  NameList<int> bondTypeNames;        
  NameList<int> angleTypeNames;       
  NameList<int> dihedralTypeNames;    
  NameList<int> improperTypeNames;    

  // Data related to quantum mechanics simulations.
  QMData *qm_data;

  // Interface to standard extra data.  These will exist and be initialized
  // to zero when init_atoms is called.

  enum dataset_flag {       
               NODATA=0x0000,
            INSERTION=0x0001, 
            OCCUPANCY=0x0002, 
              BFACTOR=0x0004, 
                 MASS=0x0008, 
               CHARGE=0x0010, 
               RADIUS=0x0020, 
               ALTLOC=0x0040, 
         ATOMICNUMBER=0x0080, 
                BONDS=0x0100, 
           BONDORDERS=0x0200,
               ANGLES=0x0400,
               CTERMS=0x0800, 
            BONDTYPES=0x1000,
           ANGLETYPES=0x2000,
            ATOMFLAGS=0x4000 };
  
  int datasetflags;
  void set_dataset_flag(int flag) { 
    datasetflags |= flag;
  }
  void unset_dataset_flag(int flag) {
    datasetflags &= (~flag);
  }
  int test_dataset_flag(int flag) {
    return (datasetflags & flag) != 0;
  }

  float *radius() { return extraflt.data("radius"); }
  float *mass() { return extraflt.data("mass");   }
  float *charge() { return extraflt.data("charge"); }
  float *beta() { return extraflt.data("beta");   }
  float *occupancy() { return extraflt.data("occupancy"); }

  // per-atom flags array
  unsigned char *flags() { return extraflg.data("flags"); }

  int radii_minmax_need_update;
  float radii_min;
  float radii_max;
  void set_radii_changed() { radii_minmax_need_update = 1; }

  void get_radii_minmax(float &min, float &max) {
    const float *r = NULL;
    if (radii_minmax_need_update && nAtoms > 0 &&
        (r = extraflt.data("radius")) != NULL) {

#if 1
      // use fast 16-byte-aligned min/max routine
      minmax_1fv_aligned(r, nAtoms, &radii_min, &radii_max);
#else
      // scalar min/max loop
      radii_min = r[0];
      radii_max = r[0];
      for (int i=0; i<nAtoms; i++) {
        if (r[i] < radii_min) radii_min = r[i];
        if (r[i] > radii_max) radii_max = r[i];
      }
#endif
      radii_minmax_need_update = 0;
    }

    min = radii_min;  
    max = radii_max;  
  }

  float *bondorders() { return extraflt.data("bondorders"); }
  void setbondorder(int atom, int bond, float order);
  float getbondorder(int atom, int bond);

  int *bondtypes() { return extraint.data("bondtypes"); }

  void setbondtype(int atom, int bond, int type);

  int  getbondtype(int atom, int bond);

  int has_structure() const { return cur_atom > 0; }

  void clear_bonds(void);

  int count_bonds(void);

#ifdef VMDWITHCARBS

  void find_small_rings_and_links(int maxpathlength, int maxringsize);
#endif

private:
  const int ID;          

  //
  // molecular structure data:
  //
  int cur_atom;          
  MolAtom *atomList;     
  int lastbonderratomid; 
  int bonderrorcount;    
 
  //
  // routines to determine components of molecular structure
  //
 
  int find_backbone(void);

  // find the residues in the molecule; return number found.
  // I look for atoms with the same resid connected to 
  // backbone atoms with the same resid (find enough backbone
  // atoms, then find atoms connected to them)
  int find_connected_backbone(IntStackHandle, int, int, int, int, int *);
  void clean_up_connection(IntStackHandle, int, int, int *);
  void find_connected_atoms_in_resid(IntStackHandle, int, int, 
     int, int, int *);
  void find_and_mark(IntStackHandle, int, int, int, int *, int *);
  int make_uniq_resids(int *flgs); 

  int find_residues(void);
  
  // This should take place after find_residues to keep
  // from mistaking a protein resname as a water resname, maybe
  void find_connected_waters(int i, char *tmp);
  int find_connected_waters2(void);

  int find_waters(void);

  void find_connected_residues(int num_residues);
  
  int find_segments(void) { return segNames.num(); }
  
  int find_connected_fragments();

  int find_fragments(void);

  void find_subfragments_cyclic(ResizeArray<Fragment *> *subfragList, int restype);  
  void find_cyclic_subfragments(ResizeArray<Fragment *> *subfragList, ResizeArray<int> *subfragCyclic);
 
  void find_connected_subfragment(int resnum, int fragnum, char *flgs, 
         int endatom, int altendatom, int alt2endatom, int alt3endatom,
         int restype, 
         ResizeArray<Fragment *> *subfragList);

  void find_subfragments(int startatom, int altstartatom, int alt2startatom,
    int endatom, int altendatom, int alt2endatom, int alt3endatom,
    int restype, ResizeArray<Fragment *> *subfragList);

  void find_subfragments_topologically(int restype, ResizeArray<Fragment *> *subfragList, int endatom, int altendatom, int alt2endatom, int alt3endatom);

#if defined(VMDFASTRIBBONS)
  void calculate_ribbon_controlpoints();
#endif

#ifdef VMDWITHCARBS

   int find_small_rings(int maxringsize);
   int find_back_edges(ResizeArray<int> &back_edge_src, ResizeArray<int> &back_edge_dest);
   int find_connected_subgraph_back_edges(int atomid, ResizeArray<int> &back_edge_src, ResizeArray<int> &back_edge_dest,
                                          int *intree_parents);
   int find_small_rings_from_back_edges(int maxringsize, ResizeArray<int> &back_edge_src, ResizeArray<int> &back_edge_dest);
   int find_small_rings_from_partial(SmallRing *ring, int maxringsize, inthash_t *used_edges, inthash_t *used_atoms);
   int get_edge_key(int edge_src, int edge_dest);
   
   // orientate small rings
   void orientate_small_rings(int maxringsize);
   void orientate_small_ring(SmallRing &ring,int maxringsize);
   
   // find links between small rings
   int find_orientated_small_ring_linkages(int maxpathlength,int maxringsize);
   int find_linkages_for_ring_from_partial(LinkagePath &lp, int maxpathlength, inthash_t *atom_to_ring, inthash_t *multi_ring_atoms, inthash_t *used_atoms);
#endif


protected:
  char *moleculename;  
  int need_find_bonds; 

public:
  // constructor; just sets most things to null.  Derived classes must put
  // in structure in 'create' routine.  Typical sequence of creating a
  // molecule should be:
  //    mol = new Molecule(....)
  //    ( specify parameters for creation )
  //    mol->create();  ... return success
  //    mol->analyze(); ... find information about the structure
  BaseMolecule(int);      
  virtual ~BaseMolecule(void); 

  //
  // routines to develop molecular structure
  //
  
  int init_atoms(int n); 

  void find_bonds_from_timestep() { need_find_bonds = 1; }
  void find_unique_bonds_from_timestep() { need_find_bonds = 2; }

  int add_atoms(int natoms, 
                const char *name, const char *type, int atomicnumber, 
                const char *resname, int resid,
                const char *chainid, const char *segname,
                const char *insertion = (char *) " ", 
                const char *altloc = "");

  int add_bond(int, int, float, int, int = ATOMNORMAL);

  int add_bond_dupcheck(int, int, float, int);

  void clear_angles(void) { angles.clear(); angleTypes.clear(); }

  int num_angles() { return int(angles.num() / 3); }

  int add_angle(int a, int b, int c, int type=-1);

  int set_angletype(int a, int type);

  int get_angletype(int a);

  void clear_dihedrals(void) { dihedrals.clear(); dihedralTypes.clear(); }

  int num_dihedrals() { return int(dihedrals.num() / 4); }

  int add_dihedral(int a, int b, int c, int d, int type=-1);
    
  int set_dihedraltype(int d, int type);

  int get_dihedraltype(int d);

  void clear_impropers(void) { impropers.clear(); improperTypes.clear(); }

  int num_impropers() { return int(impropers.num() / 4); }

  int add_improper(int a, int b, int c, int d, int type=-1);
  
  int set_impropertype(int i, int type);
  int get_impropertype(int i);

  void clear_cterms() {cterms.clear();}

  int num_cterms() { return int(cterms.num() / 8); }

  void add_cterm(int a, int b, int c, int d, 
                 int e, int f, int g, int h)  {
    cterms.append(a); cterms.append(b); cterms.append(c); cterms.append(d); 
    cterms.append(e); cterms.append(f); cterms.append(g); cterms.append(h); 
  }
  
  // (By this time, the molecule is on the MoleculeList!)
  void analyze(void);

  //
  // interfaces to add/query/delete instance transformation matrices
  //
  void add_instance(Matrix4 & inst) { instances.append(inst); }
  int num_instances(void) { return int(instances.num()); }
  void clear_instances(void) { instances.clear(); }


  //
  // query info about the molecule
  //
  int id(void) const { return ID; } 
  const char *molname() const {return moleculename; } 

  // Return the Nth atom, residue, and fragment.  All assume correct index 
  // and that the structure has been initialized (for speed).
  MolAtom *atom(int n) { return atomList+n; } 
  Residue *residue(int);                      
  Fragment *fragment(int);                    

  // return the residue or fragment in which the given atom is located.
  Residue *atom_residue(int);       
  Fragment *atom_fragment(int);     


  int find_atom_in_residue(int atomnameindex, int residue) {
    const ResizeArray<int> &atoms = residueList[residue]->atoms;
    int num = int(atoms.num());
    for (int i=0; i<num; i++) {
      if (atom(atoms[i])->nameindex == atomnameindex) return atoms[i];
    }
    return -3;
  }

  int find_atom_in_residue(const char *atomname, int residue);


  float default_charge(const char *);
  float default_mass(const char *);
  float default_radius(const char *);
  float default_occup(void) { return 1.0; }
  float default_beta(void) { return 0.0; }

  void add_volume_data(const char *name, const float *o, 
    const float *xa, const float *ya, const float *za, int x, int y, int z,
    float *voldata, float *gradient=NULL, float *variance=NULL); 

  void add_volume_data(const char *name, const double *o, 
    const double *xa, const double *ya, const double *za, int x, int y, int z,
    float *voldata); 

  void remove_volume_data(int idx);

  int num_volume_data(); 
  const VolumetricData *get_volume_data(int); 
  VolumetricData *modify_volume_data(int); 
  void compute_volume_gradient(VolumetricData *);  

protected:
  ResizeArray<VolumetricData *>volumeList;    
};

// Hydrogen atom name detection macro
#define IS_HYDROGEN(s) (s[0] == 'H' || (isdigit(s[0]) && s[1] == 'H' ))

#endif

---