psi::extrema::zmat Class Reference

Derived class for z-matrix coordinates. More...

#include <zmat.h>

Inheritance diagram for psi::extrema::zmat:

Public Member Functions
	zmat ()
	zmat constructor
	~zmat ()
	zmat destructor
void	optimize ()
	z-matrix optimization driver
Protected Member Functions
void	compute_B (void)
	Computes the B matrix.
void	cart_to_internal (double **)
void	print_internals (void)
	Prints z-matrix.
void	initial_Hi (void)
	Froms initial guess for inverse hessian.
void	grad_trans ()
	Transforms gradients from cartesians to internals.
void	parse_input ()
	Performs input parsing for z-matrix coordinates.
void	write_chkpt ()
	Writes geometries to chkpt.
void	newton_step ()
	Computes newton-raphson z-matrix optimization step.
void	back_transform ()
	Computes cartesians corresponding to the current z-matrix.
void	print_carts (double)
	Prints cartesians.
void	print_c_grads ()
	Prints cartesian gradients (Hartree/Bohr).
Protected Attributes
struct z_entry *	z_geom
simple *	simples
int *	first_unique
int *	pos_neg_pairs
double *	fcoord_old
double	bond_lim
double	angle_lim

---

Detailed Description

Derived class for z-matrix coordinates.

A top level class which allows z-matrix manipulations by driving member functions of the abstract classes. Currently optimization is the only available method.

Definition at line 22 of file zmat.h.

---

Constructor & Destructor Documentation

zmat::zmat

(

)

zmat constructor

Constructor for the top-level z-matrix derived class. Performs setup of z-matrix coordinates which includes determining optimized coordinates and positive/negative torsion pairs. The internals dummy constructor is called initially and the actual constructor named constructor_internals is called after the number of optimized coordinates is determined.

Definition at line 29 of file bin/extrema/Zmat/zmat.cc.

References psi::extrema::coord_base::angle_abort, ANGLE_TYPE, psi::extrema::coord_base_carts::c_grads, psi::extrema::coord_base_carts::carts, chkpt_rd_felement(), chkpt_rd_fgeom(), chkpt_rd_nallatom(), chkpt_rd_zmat(), psi::extrema::coord_base::coords, psi::extrema::internals::fcoords, psi::extrema::internals::fcoords_old, psi::extrema::coord_base::felement, first_unique, psi::extrema::internals::fnum_coords, free_matrix(), psi::extrema::simple::get_val(), init_array(), init_int_array(), psi::extrema::internals::mem_alloc(), NEAR_180, NOT_180, psi::extrema::coord_base_carts::num_atoms, psi::extrema::coord_base::num_coords, psi::extrema::coord_base_carts::num_entries, parse_input(), pos_neg_pairs, POS_NEG_TORS, psi::extrema::simple::set_label(), psi::extrema::simple::set_simple(), simples, TORS_TYPE, and z_geom.

             : internals()
   {

  int i, j, pos, dummy;        

  zmat::parse_input();

  /*read z_mat and cartesians from chkpt*/
  num_entries = chkpt_rd_nallatom();
  z_geom = chkpt_rd_zmat();
  char **temp_felement;
  felement = chkpt_rd_felement();

  dummy=0;
  for(i=0;i<num_entries;++i) 
      if(!strncmp(felement[i],"X\0",2) )
          dummy=1;
  if(dummy) {

      free(carts);
      double** cart_temp;
      carts = init_array(3*num_entries);
      cart_temp = chkpt_rd_fgeom();

      for(i=0;i<(3*num_entries);++i) 
          carts[i] = cart_temp[0][i];
      free_matrix(cart_temp,1);

      double* ctemp;
      ctemp = init_array(3*num_atoms);
      for(i=0;i<3*num_atoms;++i)
              ctemp[i] = c_grads[i];
      free(c_grads);
      c_grads = init_array(3*num_entries);
      pos=0;
      for(i=0;i<num_entries;++i) {
          if(strncmp(felement[i],"X\0",2) ) {
              c_grads[3*i] = ctemp[3*pos];
              c_grads[3*i+1] = ctemp[3*pos+1];
              c_grads[3*i+2] = ctemp[3*pos+2];
              ++pos;
          }
      }  
      free(ctemp);
  }

  switch(num_entries) {
    case 2: fnum_coords = 1; break;
    case 3: fnum_coords = 3; break;
    default: fnum_coords =  (num_entries*3-6); break;
  } 

  /*write z_mat to the array of simple_internal*/
  simples = (simple*) malloc(fnum_coords*sizeof(simple));  
  for(i=1;i<num_entries;++i) {
      if( i==1 ) {
          simples[0].set_simple(0,z_geom[1].bond_val,
                                2,z_geom[1].bond_atom,-1,-1,
                                z_geom[1].bond_opt);
      }
      else if( i==2 ) {
          simples[1].set_simple(0,z_geom[2].bond_val,
                                3,z_geom[2].bond_atom,-1,-1,
                                z_geom[2].bond_opt);
          simples[2].set_simple(1,z_geom[2].angle_val *_pi/180.0,
                                3,z_geom[2].bond_atom,z_geom[2].angle_atom,-1,
                                z_geom[2].angle_opt);
          j=3;
        }
      else if( i>2 ) {
          simples[j].set_simple(0,z_geom[i].bond_val,i+1,z_geom[i].bond_atom,
                                -1,-1,z_geom[i].bond_opt);
          simples[j+1].set_simple(1,z_geom[i].angle_val*_pi/180.0,
                                  i+1,z_geom[i].bond_atom,
                                  z_geom[i].angle_atom,-1,z_geom[i].angle_opt);
          simples[j+2].set_simple(2,z_geom[i].tors_val*_pi/180.0,
                                  i+1,z_geom[i].bond_atom,
                                  z_geom[i].angle_atom,z_geom[i].tors_atom,
                                  z_geom[i].tors_opt);
          j+=3;
        }
  }

  int p=0;
  for(i=1;i<num_entries;++i) {
      if(i==1) {
          simples[p].set_label(z_geom[1].bond_label);
          ++p;
      }
      else if( i==2) {
          simples[p].set_label(z_geom[2].bond_label);
          simples[p+1].set_label(z_geom[2].angle_label );
          p+=2;
      }
      else if (i>2) {
          simples[p].set_label(z_geom[i].bond_label);
          simples[p+1].set_label(z_geom[i].angle_label);
          simples[p+2].set_label(z_geom[i].tors_label);
          p+=3;
      }
  }
  
  /*find first instance of each unique coordinate*/
  int there;
  first_unique = (int *) malloc(fnum_coords*sizeof(int));
  for(i=0;i<fnum_coords;++i) {
      there=0;
      for(j=0;j<i;++j) 
          if( !strcmp( simples[i].get_label(), simples[j].get_label() )) {
              first_unique[i]=0;
              ++there;
          }
      if( (there==0) || !strcmp(simples[i].get_label(),"") )
          first_unique[i]=1;      
  }

  num_coords = 0;
  for(i=0;i<fnum_coords;++i) {
      if(simples[i].get_opt() && first_unique[i]) {
          ++num_coords;
      }
  }
  if(!num_coords)
      punt("No coordinates to optimize");

  /*allocate memory now that optimized coordinate number is known*/
  internals::mem_alloc();

  for(i=0;i<fnum_coords;++i)
      fcoords[i] = fcoords_old[i] = simples[i].get_val();

  p=0;
  for(i=0;i<fnum_coords;++i) 
      if(simples[i].get_opt() && first_unique[i]) {  
          coords[p] = simples[i].get_val();
          ++p;
      }
          
  /*find positive/negative torsion pairs*/
  int *is_torsion;
  is_torsion = init_int_array(num_coords);
  p=0;
  for(i=0;i<fnum_coords;++i) {
      if(simples[i].get_opt() && first_unique[i]) {
          if(simples[i].get_type()==TORS_TYPE) 
              is_torsion[p] = 1;
          else is_torsion[p] = 0;
          ++p;
      }
  }

  pos_neg_pairs = (int *) malloc(num_coords*sizeof(int));
  for(i=0;i<num_coords;++i)
      pos_neg_pairs[i] = 0;
  p=1;
  for(i=0;i<num_coords;++i) {
      if((pos_neg_pairs[i]==0) && is_torsion[i]) 
          for(j=(i+1);j<num_coords;++j) 
              if((pos_neg_pairs[j]==0) && is_torsion[j]) 
                  if( fabs(coords[i]+coords[j]) < POS_NEG_TORS ) {
                      pos_neg_pairs[i] = pos_neg_pairs[j] = p;
                      ++p;
                  }
  }
  
  free(is_torsion);
        
  /* test angles for extreme values, abort if hopeless */
  if(angle_abort)
      for(i=0;i<fnum_coords;++i) {
          if( simples[i].get_type() == ANGLE_TYPE ) {
              if( fabs(simples[i].get_val()) > NEAR_180*_pi/180.0 ) 
                  punt("Simple valence angle near 180 degrees");
          }
          else if (simples[i].get_type() == TORS_TYPE) 
              if( (fabs(simples[i].get_val()) > NEAR_180*_pi/180.0 ) &&
                  (fabs(simples[i].get_val()) < NOT_180*_pi/180.0) )
                  punt("Simple torsion near 180 degrees");
      }
  
  return;
}

psi::extrema::zmat::~zmat

(

)

[inline]

zmat destructor

Frees allocated memory.

Definition at line 68 of file zmat.h.

References first_unique, pos_neg_pairs, simples, and z_geom.

            {
        int i; 
        free(z_geom); 
        free(simples);
        free(first_unique);
        free(pos_neg_pairs);
        return; }

---

Member Function Documentation

void zmat::compute_B

(

void

)

[protected, virtual]

Computes the B matrix.

Interface for internals::B_row...() functions.

Implements psi::extrema::internals.

Definition at line 25 of file zmat_transform.cc.

References psi::extrema::internals::B, psi::extrema::internals::B_red, psi::extrema::internals::B_row_angle(), psi::extrema::internals::B_row_bond(), psi::extrema::internals::B_row_tors(), psi::extrema::coord_base_carts::carts, psi::extrema::coord_base::felement, first_unique, psi::extrema::internals::fnum_coords, psi::extrema::coord_base_carts::masses, psi::extrema::coord_base_carts::num_entries, simples, and psi::extrema::internals::u.

Referenced by optimize().

                     {  

    int i, j, pos=0;
    double *B_row0, *B_row1, *B_row2;

  for(i=1;i<num_entries;++i) {
      if(i==1) {
          B_row0 = B_row_bond(carts, i, simples[pos].get_bond()-1);
          B[pos] = B_row0;
          ++pos;
        }
      if(i==2) {
          B_row0 = B_row_bond(carts, i, simples[pos].get_bond()-1);
          B_row1 = B_row_angle(carts, i, simples[pos+1].get_bond()-1, 
                               simples[pos+1].get_angle()-1);
          B[pos] = B_row0;
          B[pos+1] = B_row1;
          pos += 2;
        }
      if(i>2) {
          B_row0 = B_row_bond(carts, i, simples[pos].get_bond()-1);
          B_row1 = B_row_angle(carts, i, simples[pos+1].get_bond()-1, 
                               simples[pos+1].get_angle()-1);
          B_row2 = B_row_tors(carts, i, simples[pos+2].get_bond()-1, 
                              simples[pos+2].get_angle()-1, 
                              simples[pos+2].get_tors()-1);
          B[pos] = B_row0;
          B[pos+1] = B_row1;
          B[pos+2] = B_row2;
          pos += 3;
        }
    }

  /*form u*/
  int k=0;
  for(j=0;j<num_entries;++j) {
      if(strncmp(felement[j],"X\0",2)) {
         u[3*j][3*j] = 1.0 / masses[k]; 
         u[3*j+1][3*j+1] = 1.0 /masses[k];
         u[3*j+2][3*j+2] = 1.0 / masses[k];
         ++k;
       }
       else if (!strncmp(felement[j],"X\0",2)) {
          u[3*j][3*j] = u[3*j+1][3*j+1] = u[3*j+2][3*j+2]= 1.0;
          }
  }

  /*form B_red, the reduced dimension B matrix*/

  pos=0;
  for(i=0;i<fnum_coords;++i) 
      if(first_unique[i] && simples[i].get_opt() ) {
          for(j=0;j<(3*num_entries);++j)
              B_red[pos][j] = B[i][j];
          ++pos;
      }  
      
  return;
}

void zmat::grad_trans

(

)

[protected, virtual]

Transforms gradients from cartesians to internals.

Interface for internals::grad_trans(). Gradients are calculated for redundant coordinates and their equivalency is checked.

Reimplemented from psi::extrema::internals.

Definition at line 95 of file zmat_transform.cc.

References EQUIV_GRAD, psi::extrema::internals::fgrads, first_unique, psi::extrema::internals::fnum_coords, psi::extrema::internals::grad_trans(), psi::extrema::coord_base::grads, NORMAL_PRINT, psi::extrema::coord_base::print_lvl, and simples.

Referenced by optimize().

                      {
    
    internals::grad_trans();

    int i, j, div, p=0;

    if(print_lvl>NORMAL_PRINT) {
        fprintf(outfile,"\n  Internal coordinate gradients (a.u):");
        for(i=0;i<fnum_coords;++i)
            fprintf(outfile,"\n  %i %8s: %15.10lf",
                    i+1,simples[i].get_label(),fgrads[i]);
    }
    fprintf(outfile,"\n");

    double sum;
    for(i=0;i<fnum_coords;++i) 
        if(first_unique[i] && simples[i].get_opt()) {
            sum = fgrads[i];
            div=1;
            for(j=(i+1);j<fnum_coords;++j)
                if(!strcmp(simples[i].get_label(),simples[j].get_label()) &&
                    strcmp(simples[i].get_label(),"")) {
                    sum += fgrads[j];
                    ++div;
                    if( fabs(fgrads[i]-fgrads[j]) > EQUIV_GRAD ) {
                        fprintf(outfile,"\n  WARNING: gradients %d and %d",
                                i+1,j+1);
                        fprintf(outfile," should be equal, differ by %lf",
                                fabs(fgrads[i]-fgrads[j]));
                    }
                }
            grads[p] = sum/((double) div);
            ++p;
        }
    
    return;
}

void zmat::newton_step

(

)

[protected]

Computes newton-raphson z-matrix optimization step.

Interface for math_tools::newton_step().

Definition at line 405 of file bin/extrema/Zmat/zmat.cc.

References psi::extrema::coord_base::angle_abort, angle_lim, ANGLE_TYPE, bond_lim, psi::extrema::coord_base::coord_write, psi::extrema::coord_base::coords, fcoord_old, first_unique, psi::extrema::internals::fnum_coords, psi::extrema::coord_base::grads, psi::extrema::coord_base::Hi, init_array(), NEAR_180, psi::extrema::math_tools::newton_step(), NOT_180, psi::extrema::coord_base::num_coords, pos_neg_pairs, POS_NEG_TORS, psi::extrema::coord_base::print_lvl, psi::extrema::simple::set_val(), simples, and TORS_TYPE.

Referenced by optimize().

                         {
    
    int i, j, k, p;
    double *s, num, con;

    fprintf(outfile,"\n\n  --------------------------------------");
    fprintf(outfile,"--------------------------------------\n");
    fprintf(outfile,"  Computing newton-raphson optimization step\n");
    fprintf(outfile,"  --------------------------------------");
    fprintf(outfile,"--------------------------------------\n");
    s = math_tools::newton_step(num_coords,Hi,grads);

    /* ensure pos/neg torsion angle pairs match */
    for(i=0;i<num_coords;++i) 
        if(pos_neg_pairs[i]) 
            for(j=(i+1);j<num_coords;++j) 
                if(pos_neg_pairs[i]==pos_neg_pairs[j]) {
                    
                    num = (fabs(s[i]) + fabs(s[j]))/2;
                    if(s[i] < 0)
                        s[i] = -1.0*num;
                    else s[i] = num;
                    if(s[j] < 0)
                        s[j] = -1.0*num;
                    else s[j] = num;
                    
                    if( (s[i]+s[j]) > POS_NEG_TORS ) {
                        fprintf(outfile,"\n  WARNING: positive/negative pair");
                        fprintf(outfile," %d/%d:\n",i,j);
                        fprintf(outfile,"  displacements differ by more than");
                        fprintf(outfile,"%lf radians\n", POS_NEG_TORS);
                    }
                }
    
    /* print before limit enforcement */
    if(print_lvl>1) {
        p=0;
        fprintf(outfile,"\n  Displacements before limit enforcement");
        fprintf(outfile," (angstroms and degrees):\n");
        for(i=0;i<fnum_coords;++i)
            if(simples[i].get_opt() && first_unique[i]) {
                if(simples[i].get_type()==0)
                    con = _bohr2angstroms; 
                else con = 180.0/_pi;  
                fprintf(outfile,"  %i %8s: %lf\n",
                        i+1,simples[i].get_label(),s[p]*con);
                ++p;
            }
    }
    
    /* enforce bond and angle limits */
    p=0;
    for(i=0;i<fnum_coords;++i) {
        if((simples[i].get_type()==0) && first_unique[i] && simples[i].get_opt() ) {
            if(s[p] > bond_lim) 
                s[p] = bond_lim; 
            else if (s[p] < (-1.0*bond_lim) ) 
                s[p] = (-1.0*bond_lim); 
            ++p;
        }
        else if(first_unique[i] && simples[i].get_opt() ) {
            if(s[p] > angle_lim) 
                s[p] = angle_lim; 
            else if (s[p] < (-1.0*angle_lim) ) 
                s[p] = (-1.0*angle_lim);
            ++p;
        }
    }

    /*save old coordinates to write to opt.dat and for back transform*/
    fcoord_old = init_array(fnum_coords);
    for(i=0;i<num_coords;++i) 
        coord_write[i] = coords[i];
    for(i=0;i<fnum_coords;++i)
        fcoord_old[i] = simples[i].get_val();


    /* take the step */
    for(i=0;i<num_coords;++i) 
        coords[i] += s[i];

    /*update simples*/
    p=0;
    for(i=0;i<fnum_coords;++i) {
        if(first_unique[i] && simples[i].get_opt()) {
            simples[i].set_val(coords[p]);
            ++p;
            for(j=(i+1);j<fnum_coords;++j) 
                if(!strcmp(simples[i].get_label(),simples[j].get_label()) &&
                    strcmp(simples[i].get_label(),""))
                    simples[j].set_val(simples[i].get_val());
        }
    }

    /* test for extreme angle cases */
    if(angle_abort)
        for(i=0;i<fnum_coords;++i) {
            if( simples[i].get_type() == ANGLE_TYPE ) {
                if( fabs(simples[i].get_val()) > NEAR_180*_pi/180.0 )
                    punt("Simple valence angle near 180 degrees");
            }
            else if (simples[i].get_type() == TORS_TYPE)
                if( (fabs(simples[i].get_val()) > NEAR_180*_pi/180.0 ) &&
                    (fabs(simples[i].get_val()) < NOT_180*_pi/180.0) )
                    punt("Simple torsion near 180 degrees");
        }
    
    int entry;
    fprintf(outfile,"\n  Optimization Step (angstroms and degrees):\n");
    fprintf(outfile,"\n  label    initial value   gradient (a.u.)");      
    fprintf(outfile," displacement    new value");
    fprintf(outfile,"\n  -------- --------------- ---------------"); 
    fprintf(outfile," --------------- ---------------");
    entry=0;
    p=0;
    for(i=0;i<fnum_coords;++i) {
        con = _bohr2angstroms;
        if( i==2 )
            con = 180.0/_pi;
        if( i>2 ) {
            if(entry==0)
                ++entry;
            else if(entry==1) {
                con = 180.0/_pi;
                ++entry;
            }
            else if(entry==2) {
                con = 180.0/_pi;
                entry = 0;
            }
        }
        if(simples[i].get_opt() && first_unique[i]) {
            fprintf(outfile,"\n  %-8s %15.10lf %15.10lf %15.10lf %15.10lf",
                    simples[i].get_label(),coord_write[p]*con, grads[p], 
                    s[p]*con, coords[p]*con);
            ++p;
        }
    }

    free(s);
    return;
}

void zmat::back_transform

(

)

[protected]

Computes cartesians corresponding to the current z-matrix.

Interface for internals::back_transform().

Definition at line 142 of file zmat_transform.cc.

References psi::extrema::internals::back_transform(), fcoord_old, psi::extrema::internals::fnum_coords, psi::extrema::simple::get_val(), init_array(), and simples.

Referenced by optimize().

                          {
    
    /*create full coordinate vectors*/
    double *fcoord_new;
    fcoord_new = init_array(fnum_coords);
    
    int i;
    for(i=0;i<fnum_coords;++i) { 
        fcoord_new[i] = simples[i].get_val();
    }

    internals::back_transform(fcoord_new, fcoord_old);

    free(fcoord_old);
    free(fcoord_new);
    return;
}

void zmat::print_carts

(

double

conv

)

[protected, virtual]

Prints cartesians.

Parameters:

conv

conversion factor; either 1.0 for bohr or _bohr2angstroms for angstroms

Reimplemented from psi::extrema::coord_base_carts.

Definition at line 558 of file bin/extrema/Zmat/zmat.cc.

References psi::extrema::coord_base_carts::carts, psi::extrema::coord_base::felement, free_matrix(), init_matrix(), and psi::extrema::coord_base_carts::num_entries.

Referenced by optimize().

                                    {
 
    int i, j;
    double **temp;

    temp = init_matrix(num_entries,3);
    for(i=0;i<num_entries;++i) 
        for(j=0;j<3;++j) 
            temp[i][j] = carts[3*i+j]*conv;
    if(conv==1.0)
        fprintf(outfile,"\n  Cartesian Coordinates (bohr):\n");
    else
        fprintf(outfile,"\n  Cartesian Coordinates (angstroms):\n");
     fprintf(outfile,"                       x              y         ");
    fprintf(outfile,"       z\n");
    fprintf(outfile,"                --------------- --------------- ");
    fprintf(outfile,"---------------\n");
    for(i=0;i<num_entries;++i)
        fprintf(outfile,"  %12s  %15.10lf %15.10lf %15.10lf\n",
                felement[i], temp[i][0], temp[i][1], temp[i][2]);
    free_matrix(temp,num_entries);

    return;
}

---

Field Documentation

struct z_entry* psi::extrema::zmat::z_geom [read, protected]

an array of z_entry structs, defined in chkpt.h, a z-matrix is read/written to chkpt as an array of z_entry structs

Definition at line 26 of file zmat.h.

Referenced by write_chkpt(), zmat(), and ~zmat().

simple* psi::extrema::zmat::simples [protected]

an array of simple classes, the simple concrete class provides interfaces to simple internal coordinate information

Definition at line 32 of file zmat.h.

Referenced by back_transform(), compute_B(), grad_trans(), initial_Hi(), newton_step(), print_internals(), write_chkpt(), zmat(), and ~zmat().

int* psi::extrema::zmat::first_unique [protected]

1 if a simple intenal from the full set of simples is the first unique coordinate, 0 otherwise, only unique coordinates need to be optimized

Definition at line 36 of file zmat.h.

Referenced by compute_B(), grad_trans(), initial_Hi(), newton_step(), zmat(), and ~zmat().

int* psi::extrema::zmat::pos_neg_pairs [protected]

positive negative torsion angle pairs are indicated by matching positive integers, for keeping torsion angle pairs matching

Definition at line 41 of file zmat.h.

Referenced by newton_step(), zmat(), and ~zmat().

double* psi::extrema::zmat::fcoord_old [protected]

full coordinates prior to optimization step, needed for back transformation to cartesians

Definition at line 46 of file zmat.h.

Referenced by back_transform(), and newton_step().

double psi::extrema::zmat::bond_lim [protected]

maximum bond displacement in bohr

Definition at line 49 of file zmat.h.

Referenced by newton_step(), and parse_input().

double psi::extrema::zmat::angle_lim [protected]

maximum angle displacement in radians

Definition at line 50 of file zmat.h.

Referenced by newton_step(), and parse_input().

---

The documentation for this class was generated from the following files:

• zmat.h
• bin/extrema/Zmat/zmat.cc
• zmat_transform.cc

---