xyz2pov/src/pov_writer.c

#define _FILE_OFFSET_BITS 64

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "pov_writer.h"
#include "atom_parser.h"

struct bond{
  int i;
  int j;
};

struct linked_bond_list{
  struct bond the_bond;
  struct linked_bond_list *next;
};

void make_bonds(struct coords *, int);

struct linked_bond_list *bl_head;
  
void clean_bonds(void);

void initBondList(void){
  bl_head = NULL;
}

void pov_write(FILE *out_file, struct coords *the_coords, int n_atoms,
               int d_hydrogens, int d_bonds, int d_atoms, int d_vectors){

  int i,j; /*loop counters */
  int skip_atom, skip_bond, test1, test2; /*booleans */
  int gb_atom, gbdp_atom;
  double dx, dy, dz; /* used in making the bonds */
  
  struct linked_bond_list *current_bond; /*keeps track of the linked list*/

  for(i = 0; i < n_atoms; i++){
    update_types(the_coords[i].name);
  }  
  
  if(d_atoms){
    
    fprintf(out_file,
            "//************************************************************\n"
            "// The list of atoms\n"
            "//************************************************************\n"
            "\n"
            "\n");
    
    for(i = 0; i < n_atoms; i++){
      
      skip_atom = 0;     
      
      if(!d_hydrogens){
        skip_atom = !strcmp("H", the_coords[i].name);
      }
      
      if(!skip_atom){          
        
        gb_atom = !strcmp("GB", the_coords[i].name);
        gbdp_atom = !strcmp("GBDP", the_coords[i].name);
        
        if (gb_atom) {
          fprintf(out_file, 
                  "make_%s_ellipse( %lf, %lf, %lf, %lf, %lf, %lf, %lf)\n",
                  the_coords[i].name,
                  the_coords[i].x,
                  the_coords[i].z,
                  the_coords[i].y,
                  the_coords[i].charge,
                  the_coords[i].ux,
                  the_coords[i].uz,
                  the_coords[i].uy);
        } else {
          if (gbdp_atom) {
            fprintf(out_file, 
                    "make_%s_shaded_ellipse( %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf)\n",
                    the_coords[i].name,
                    the_coords[i].x,
                    the_coords[i].z,
                    the_coords[i].y,
                    the_coords[i].charge,
                    the_coords[i].ux,
                    the_coords[i].uz,
                    the_coords[i].uy,
                    the_coords[i].vx,
                    the_coords[i].vz,
                    the_coords[i].vy);
          } else {
            fprintf(out_file, 
                    "make_%s_atom( %lf, %lf, %lf )\n",
                    the_coords[i].name,
                    the_coords[i].x,
                    the_coords[i].z,
                    the_coords[i].y);
          }
        }
      }
    }
    
    fprintf(out_file,
            "\n"
            "\n");
  }

      
  if (d_vectors) {

    fprintf(out_file,
            "//************************************************************\n"
            "// The list of vectors\n"
            "//************************************************************\n"
            "\n"
            "\n");
    
    for(i = 0; i < n_atoms; i++){
      
      if (the_coords[i].hasVector) {
        fprintf(out_file,
                "make_%s_vector(%lf, %lf, %lf, %lf, %lf, %lf)\n",
                the_coords[i].name,
                the_coords[i].x,
                the_coords[i].z,
                the_coords[i].y,
                the_coords[i].ux,
                the_coords[i].uz,
                the_coords[i].uy);
      }
    }
    
    fprintf(out_file,
            "\n"
            "\n");
  }
      
  if(d_bonds){
    
    fprintf(out_file,
            "//************************************************************\n"
            "// The list of bonds\n"
            "//************************************************************\n"
            "\n"
            "\n");
    
    if( bl_head == NULL ) make_bonds(the_coords, n_atoms);
    
    current_bond = bl_head->next;
    
    while(current_bond != NULL){
      
      skip_bond = 0;
      
      i = current_bond->the_bond.i;
      j = current_bond->the_bond.j;
      
      if(!d_hydrogens){

        test1 = !strcmp("H", the_coords[i].name);
        test2 = !strcmp("H", the_coords[j].name);

        skip_bond = (test1 || test2);
      }

      if(!skip_bond){
        
        dx = (the_coords[j].x - the_coords[i].x) / 2.0;
        dy = (the_coords[j].y - the_coords[i].y) / 2.0;
        dz = (the_coords[j].z - the_coords[i].z) / 2.0;
        
        fprintf(out_file,
               "make_%s_bond( %lf, %lf, %lf, %lf, %lf, %lf )\n",
               the_coords[i].name,
               the_coords[i].x,
               the_coords[i].z,
               the_coords[i].y,
               (the_coords[i].x + dx),
               (the_coords[i].z + dz),
               (the_coords[i].y + dy));
        
        fprintf(out_file,
               "make_%s_bond( %lf, %lf, %lf, %lf, %lf, %lf )\n",
               the_coords[j].name,
               the_coords[j].x,
               the_coords[j].z,
               the_coords[j].y,
               (the_coords[j].x - dx),
               (the_coords[j].z - dz),
               (the_coords[j].y - dy));

        fprintf(out_file, "\n");
      }
      
      current_bond = current_bond->next;
    }
  
    if( regenerateBonds )clean_bonds();
  }
}


void make_bonds(struct coords *the_coords, int n_atoms){
  
  int i, j; /*counters */
  struct linked_bond_list *temp_bond; /*bond place holder */
  
  const double bond_fudge = 1.12; // a fudge factor
  struct atom type_i, type_j; /* holds the atom types */

  int test; /* booleans */
  double dx, dy, dz, dr2, dcv, dcv2; // used to determine bond existence 

  
  bl_head = (struct linked_bond_list *)malloc(sizeof(struct linked_bond_list));
  bl_head->next = NULL;

  for(i = 0; i < (n_atoms - 1); i++){
    
    for(j = (i+1); j < n_atoms; j++){
      
      dx = the_coords[j].x - the_coords[i].x;
      dy = the_coords[j].y - the_coords[i].y;
      dz = the_coords[j].z - the_coords[i].z;
      
      dr2 = dx * dx + dy * dy + dz * dz;
      
      test = !findAtomType(the_coords[i].name, &type_i);
      if(test){
        fprintf(stderr, "Atom Type %s, not found!\n",
                the_coords[i].name);
        exit(8);
      }

      test = !findAtomType(the_coords[j].name, &type_j);
      if(test){
        fprintf(stderr, "Atom Type %s, not found!\n",
                the_coords[j].name);
        exit(8);
      }
      
      
      dcv = bond_fudge * (type_i.covalentRadii + type_j.covalentRadii);
      dcv2 = dcv * dcv;

      if(dr2 <= dcv2){
        
        temp_bond = 
          (struct linked_bond_list *)malloc(sizeof(struct linked_bond_list));
        
        bl_head->the_bond.i = i;
        bl_head->the_bond.j = j;
        
        temp_bond->next = bl_head;
        bl_head = temp_bond;
      }
    }
  }
}
        

void clean_bonds(){
  struct linked_bond_list *current_bond;
  struct linked_bond_list *next_bond; /* place holders */
  

    current_bond = bl_head->next;

    while(current_bond != NULL){

      next_bond = current_bond->next;
      free(current_bond);
      current_bond = next_bond;
    }
    
    bl_head->next = NULL;
    free( bl_head );
    bl_head = NULL;
}


void make_header_macros(FILE *out_file){

  struct linked_atom *type_list;  // list of all atom types
  struct linked_atom *current_type; // current atom type

  char *name;
  double red, green, blue;
  double radius;
  
  type_list = get_type_list();
  current_type = type_list->next;
  
  while(current_type != NULL){
    
    name = current_type->myAtom.name;
    radius = current_type->myAtom.vanDerWallRadii;
    red = ((double)current_type->myAtom.red) / 255.0;
    green =  ((double)current_type->myAtom.green) / 255.0;
    blue =  ((double)current_type->myAtom.blue) / 255.0;

    
    fprintf(out_file,
            "//****************************************************\n"
            "// DEFINE %s MACROS\n" 
            "//****************************************************\n"
            "\n"
            "#macro make_%s_bond "
            "(end_1x, end_1y, end_1z, end_2x, end_2y, end_2z)\n"
            "\n"
            "  #local x1 = end_1x;\n"
            "  #local y1 = end_1y;\n"
            "  #local z1 = end_1z;\n"
            "  #local x2 = end_2x;\n"
            "  #local y2 = end_2y;\n"
            "  #local z2 = end_2z;\n"
            "\n"
            "  #if(ROTATE)\n"
            "    #local x1_new = rotatePointX + A11 * (x1-rotatePointX) + A12 * (y1-rotatePointY) + A13 * (z1-rotatePointZ);\n"
            "    #local y1_new = rotatePointY + A21 * (x1-rotatePointX) + A22 * (y1-rotatePointY) + A23 * (z1-rotatePointZ);\n"
            "    #local z1_new = rotatePointZ + A31 * (x1-rotatePointX) + A32 * (y1-rotatePointY) + A33 * (z1-rotatePointZ);\n"
            "\n"
            "    #local x2_new = rotatePointX + A11 * (x2-rotatePointX) + A12 * (y2-rotatePointY) + A13 * (z2-rotatePointZ);\n"
            "    #local y2_new = rotatePointY + A21 * (x2-rotatePointX) + A22 * (y2-rotatePointY) + A23 * (z2-rotatePointZ);\n"
            "    #local z2_new = rotatePointZ + A31 * (x2-rotatePointX) + A32 * (y2-rotatePointY) + A33 * (z2-rotatePointZ);\n"
            "\n"
            "  #else\n"
            "    #local x1_new = x1;"
            "    #local y1_new = y1;"
            "    #local z1_new = z1;"
            "\n"
            "    #local x2_new = x2;"
            "    #local y2_new = y2;"
            "    #local z2_new = z2;"
            "\n"
            "  #end\n"
            "\n"
            "  cylinder{\n"
            "    < x1_new, y1_new, z1_new >,\n"
            "    < x2_new, y2_new, z2_new >,\n"
            "    BOND_RADIUS\n"
            "    texture{\n"
            "      pigment{ rgb < %lf, %lf, %lf > }\n"
            "      finish{\n"
            "        ambient .2\n"
            "        diffuse .6\n"
            "        specular 1\n"
            "        roughness .001\n"
            "        metallic\n"
            "      }\n"
            "    }\n"
            "  }\n"
            "#end\n"
            "#macro make_%s_atom "
            "(center_x, center_y, center_z)\n"
                    "\n"
            "  #local x1 = center_x;\n"
            "  #local y1 = center_y;\n"
            "  #local z1 = center_z;\n"
            "\n"
            "  #if(ROTATE)\n"
            "\n"
            "    #local x1_new = rotatePointX + A11 * (x1-rotatePointX) + A12 * (y1-rotatePointY) + A13 * (z1-rotatePointZ);\n"
            "    #local y1_new = rotatePointY + A21 * (x1-rotatePointX) + A22 * (y1-rotatePointY) + A23 * (z1-rotatePointZ);\n"
            "    #local z1_new = rotatePointZ + A31 * (x1-rotatePointX) + A32 * (y1-rotatePointY) + A33 * (z1-rotatePointZ);\n"
            "\n"
            "  #else\n"
            "\n"
            "    #local x1_new = x1;"
            "    #local y1_new = y1;"
            "    #local z1_new = z1;"
            "\n"
            "  #end\n"
            "\n"
            "  sphere{\n"
            "    < x1_new, y1_new, z1_new >,\n"
            "    ATOM_SPHERE_FACTOR * %lf\n"
            "    texture{\n"
            "      pigment{ rgb < %lf, %lf, %lf > }\n"
            "      finish{\n"
            "        ambient .2\n"
            "        diffuse .6\n"
            "        specular 1\n"
            "        roughness .001\n"
            "        metallic\n"
            "      }\n"
            "    }\n"
            "  }\n"
            "#end\n"
            "#macro make_%s_ellipse "
            "(center_x, center_y, center_z, ecc, u_x, u_y, u_z)\n"
                    "\n"
            "  #local x1 = center_x;\n"
            "  #local y1 = center_y;\n"
            "  #local z1 = center_z;\n"
            "  #local x2 = u_x;\n"
            "  #local y2 = u_y;\n"
            "  #local z2 = u_z;\n"
            "\n"
            "  #if(ROTATE)\n"
            "\n"
            "    #local x1_new = rotatePointX + A11 * (x1-rotatePointX) + A12 * (y1-rotatePointY) + A13 * (z1-rotatePointZ);\n"
            "    #local y1_new = rotatePointY + A21 * (x1-rotatePointX) + A22 * (y1-rotatePointY) + A23 * (z1-rotatePointZ);\n"
            "    #local z1_new = rotatePointZ + A31 * (x1-rotatePointX) + A32 * (y1-rotatePointY) + A33 * (z1-rotatePointZ);\n"
            "\n"
            "    #local x2_new = rotatePointX + A11 * (x2-rotatePointX) + A12 * (y2-rotatePointY) + A13 * (z2-rotatePointZ);\n"
            "    #local y2_new = rotatePointY + A21 * (x2-rotatePointX) + A22 * (y2-rotatePointY) + A23 * (z2-rotatePointZ);\n"
            "    #local z2_new = rotatePointZ + A31 * (x2-rotatePointX) + A32 * (y2-rotatePointY) + A33 * (z2-rotatePointZ);\n"
            "\n"
            "  #else\n"
            "\n"
            "    #local x1_new = x1;"
            "    #local y1_new = y1;"
            "    #local z1_new = z1;"
            "\n"
            "    #local x2_new = x2;"
            "    #local y2_new = y2;"
            "    #local z2_new = z2;"
            "\n"
            "  #end\n"
            "\n"
            "    #local myUlen = sqrt(x2_new*x2_new + y2_new*y2_new + z2_new*z2_new);\n"
            "    #local uux = x2_new / myUlen;\n"
            "    #local uuy = y2_new / myUlen;\n"
            "    #local uuz = z2_new / myUlen;\n"
            "    #local myTheta = -degrees(acos(uuz));\n"
            "    #local myPsi = -degrees(atan(uux/uuy));\n"
            "    #local myScale = ATOM_SPHERE_FACTOR * %lf;\n"
            "\n"
            "  sphere{\n"
            "    < 0, 0, 0 >, 1\n"
            "    texture{\n"
            "      pigment{\n"
            "        average\n"
            "        pigment_map{\n"
            "          [1.0 grad1]\n"
            "          [1.0 grad2]\n"
            "          [1.0 grad3]\n"
            "          [5.0 gradz]\n"
            "        }\n"
            "      }\n"
            "      finish{\n"
            "        ambient .2\n"
            "        diffuse .6\n"
            "        specular 1\n"
            "        roughness .001\n"
            "        metallic\n"
            "      }\n"
            "    }\n"
            "    scale<myScale,myScale,ecc*myScale>\n"
            "    rotate<myTheta,0,myPsi>\n"
            "    translate< x1_new, y1_new, z1_new>\n"
            "  }\n"
            "#end\n"
            "#macro make_%s_shaded_ellipse "
            "(center_x, center_y, center_z, ecc, u_x, u_y, u_z, v_x, v_y, v_z)\n"
                    "\n"
            "  #local x1 = center_x;\n"
            "  #local y1 = center_y;\n"
            "  #local z1 = center_z;\n"
            "  #local x2 = u_x;\n"
            "  #local y2 = u_y;\n"
            "  #local z2 = u_z;\n"
            "  #local x3 = v_x;\n"
            "  #local y3 = v_y;\n"
            "  #local z3 = v_z;\n"
            "\n"
            "  #if(ROTATE)\n"
            "\n"
            "    #local x1_new = rotatePointX + A11 * (x1-rotatePointX) + A12 * (y1-rotatePointY) + A13 * (z1-rotatePointZ);\n"
            "    #local y1_new = rotatePointY + A21 * (x1-rotatePointX) + A22 * (y1-rotatePointY) + A23 * (z1-rotatePointZ);\n"
            "    #local z1_new = rotatePointZ + A31 * (x1-rotatePointX) + A32 * (y1-rotatePointY) + A33 * (z1-rotatePointZ);\n"
            "\n"
            "    #local x2_new = rotatePointX + A11 * (x2-rotatePointX) + A12 * (y2-rotatePointY) + A13 * (z2-rotatePointZ);\n"
            "    #local y2_new = rotatePointY + A21 * (x2-rotatePointX) + A22 * (y2-rotatePointY) + A23 * (z2-rotatePointZ);\n"
            "    #local z2_new = rotatePointZ + A31 * (x2-rotatePointX) + A32 * (y2-rotatePointY) + A33 * (z2-rotatePointZ);\n"
            "\n"
            "    #local x3_new = rotatePointX + A11 * (x3-rotatePointX) + A12 * (y3-rotatePointY) + A13 * (z3-rotatePointZ);\n"
            "    #local y3_new = rotatePointY + A21 * (x3-rotatePointX) + A22 * (y3-rotatePointY) + A23 * (z3-rotatePointZ);\n"
            "    #local z3_new = rotatePointZ + A31 * (x3-rotatePointX) + A32 * (y3-rotatePointY) + A33 * (z3-rotatePointZ);\n"
            "\n"
            "  #else\n"
            "\n"
            "    #local x1_new = x1;"
            "    #local y1_new = y1;"
            "    #local z1_new = z1;"
            "\n"
            "    #local x2_new = x2;"
            "    #local y2_new = y2;"
            "    #local z2_new = z2;"
            "\n"
            "    #local x3_new = x3;"
            "    #local y3_new = y3;"
            "    #local z3_new = z3;"
            "\n"
            "  #end\n"
            "\n"
            "    #local myUlen = sqrt(x2_new*x2_new + y2_new*y2_new + z2_new*z2_new);\n"
            "    #local uux = x2_new / myUlen;\n"
            "    #local uuy = y2_new / myUlen;\n"
            "    #local uuz = z2_new / myUlen;\n"
            "    #local myVlen = sqrt(x3_new*x3_new + y3_new*y3_new + z3_new*z3_new);\n"
            "    #local vvx = x3_new / myVlen;\n"
            "    #local vvy = y3_new / myVlen;\n"
            "    #local vvz = z3_new / myVlen;\n"
            "\n"
            "    #local myTheta = degrees(acos(uuz));\n"
            "    #local myPsi = -degrees(atan(uux/uuy));\n"
            "    #local myPhi = degrees(acos(vvz));\n"
            "    #local myScale = ATOM_SPHERE_FACTOR * %lf;\n"
            "\n"
            "  sphere{\n"
            "    < 0, 0, 0 >, 1\n"
            "    texture{\n"
            "      pigment{\n"
            "        average\n"
            "        pigment_map{\n"
            "          [1.0 grad1]\n"
            "          [1.0 grad2]\n"
            "          [1.0 grad3]\n"
            "          [5.0 gradz]\n"
            "        }\n"
            "      }\n"
            "      finish{\n"
            "        ambient .2\n"
            "        diffuse .6\n"
            "        specular 1\n"
            "        roughness .001\n"
            "        metallic\n"
            "      }\n"
            "    }\n"
            "    scale<myScale,myScale,ecc*myScale>\n"
            "    rotate<myTheta,myPhi,myPsi>\n"
            "    translate< x1_new, y1_new, z1_new>\n"
            "  }\n"
            "#end\n"
            "#macro make_%s_vector "
            "(center_x, center_y, center_z, ux, uy, uz)\n"
            "\n"
            "  #local vx = VECTOR_SCALE * ux;\n"
            "  #local vy = VECTOR_SCALE * uy;\n"
            "  #local vz = VECTOR_SCALE * uz;\n"
            "  #local x1 = center_x - 0.5 * vx;\n"
            "  #local y1 = center_y - 0.5 * vy;\n"
            "  #local z1 = center_z - 0.5 * vz;\n"
            "  #local x2 = center_x + 0.5 * vx;\n"
            "  #local y2 = center_y + 0.5 * vy;\n"
            "  #local z2 = center_z + 0.5 * vz;\n"
            "  #local v2 = vx*vx + vy*vy + vz*vz;\n"
            "  #local vl  = sqrt(v2);\n"
            "  #local x3 = x1 + vx * (1.0 - CONE_FRACTION);\n"
            "  #local y3 = y1 + vy * (1.0 - CONE_FRACTION);\n"
            "  #local z3 = z1 + vz * (1.0 - CONE_FRACTION);\n"
            "\n"
            "  #if(ROTATE)\n"
            "    #local x1_new = rotatePointX + A11 * (x1-rotatePointX) + A12 * (y1-rotatePointY) + A13 * (z1-rotatePointZ);\n"
            "    #local y1_new = rotatePointY + A21 * (x1-rotatePointX) + A22 * (y1-rotatePointY) + A23 * (z1-rotatePointZ);\n"
            "    #local z1_new = rotatePointZ + A31 * (x1-rotatePointX) + A32 * (y1-rotatePointY) + A33 * (z1-rotatePointZ);\n"
            "\n"
            "    #local x2_new = rotatePointX + A11 * (x2-rotatePointX) + A12 * (y2-rotatePointY) + A13 * (z2-rotatePointZ);\n"
            "    #local y2_new = rotatePointY + A21 * (x2-rotatePointX) + A22 * (y2-rotatePointY) + A23 * (z2-rotatePointZ);\n"
            "    #local z2_new = rotatePointZ + A31 * (x2-rotatePointX) + A32 * (y2-rotatePointY) + A33 * (z2-rotatePointZ);\n"
            "\n"
            "    #local x3_new = rotatePointX + A11 * (x3-rotatePointX) + A12 * (y3-rotatePointY) + A13 * (z3-rotatePointZ);\n"
            "    #local y3_new = rotatePointY + A21 * (x3-rotatePointX) + A22 * (y3-rotatePointY) + A23 * (z3-rotatePointZ);\n"
            "    #local z3_new = rotatePointZ + A31 * (x3-rotatePointX) + A32 * (y3-rotatePointY) + A33 * (z3-rotatePointZ);\n"
            "\n"
            "  #else\n"
            "    #local x1_new = x1;"
            "    #local y1_new = y1;"
            "    #local z1_new = z1;"
            "\n"
            "    #local x2_new = x2;"
            "    #local y2_new = y2;"
            "    #local z2_new = z2;"
            "\n"
            "    #local x3_new = x3;"
            "    #local y3_new = y3;"
            "    #local z3_new = z3;"
            "\n"
            "  #end\n"
            "\n"
            "  cylinder{\n"
            "    < x1_new, y1_new, z1_new >,\n"
            "    < x3_new, y3_new, z3_new >,\n"
            "    STICK_RADIUS\n"
            "    texture{\n"
            "      pigment{ rgb < %lf, %lf, %lf > }\n"
            "      finish{\n"
            "        ambient .2\n"
            "        diffuse .6\n"
            "        specular 1\n"
            "        roughness .001\n"
            "        metallic\n"
            "      }\n"
            "    }\n"
            "  }\n"
            "  cone{\n"
            "    < x2_new, y2_new, z2_new >, 0.0\n"
            "    < x3_new, y3_new, z3_new >, CONE_RADIUS\n"
            "    texture{\n"
            "      pigment{ rgb < %lf, %lf, %lf > }\n"
            "      finish{\n"
            "        ambient .2\n"
            "        diffuse .6\n"
            "        specular 1\n"
            "        roughness .001\n"
            "        metallic\n"
            "      }\n"
            "    }\n"
            "  }\n"
            "#end\n"
            "\n"
            "\n",
            name,
            name,
            red, green, blue,
            name,
            radius,
            red, green, blue,
            name,
            radius,
            name,
            radius,
            name,
            red, green, blue,
            red, green, blue);
    
    current_type = current_type->next;
  }
}
Revision:	2750
Committed:	Fri May 12 19:24:30 2006 UTC (18 years, 11 months ago) by tim
Content type:	text/plain
File size:	19711 byte(s)
Log Message:	adding support for Gay-Berne atom
#	User	Rev	Content
1	mmeineke	1095	#define _FILE_OFFSET_BITS 64
2
3	mmeineke	60	#include <stdio.h>
4			#include <stdlib.h>
5			#include <string.h>
6			#include <math.h>
7
8			#include "pov_writer.h"
9			#include "atom_parser.h"
10
11			struct bond{
12			int i;
13			int j;
14			};
15
16			struct linked_bond_list{
17			struct bond the_bond;
18			struct linked_bond_list *next;
19			};
20
21			void make_bonds(struct coords *, int);
22
23			struct linked_bond_list *bl_head;
24
25			void clean_bonds(void);
26
27	mmeineke	508	void initBondList(void){
28			bl_head = NULL;
29			}
30	mmeineke	60
31			void pov_write(FILE out_file, struct coords the_coords, int n_atoms,
32	gezelter	864	int d_hydrogens, int d_bonds, int d_atoms, int d_vectors){
33	mmeineke	60
34			int i,j; /loop counters /
35			int skip_atom, skip_bond, test1, test2; /booleans /
36	tim	2750	int gb_atom, gbdp_atom;
37	mmeineke	60	double dx, dy, dz; /* used in making the bonds */
38
39			struct linked_bond_list current_bond; /keeps track of the linked list*/
40
41			for(i = 0; i < n_atoms; i++){
42			update_types(the_coords[i].name);
43			}
44
45			if(d_atoms){
46
47			fprintf(out_file,
48			"//************************************************************\n"
49			"// The list of atoms\n"
50			"//************************************************************\n"
51			"\n"
52			"\n");
53	gezelter	864
54	mmeineke	60	for(i = 0; i < n_atoms; i++){
55
56	gezelter	864	skip_atom = 0;
57	mmeineke	60
58			if(!d_hydrogens){
59			skip_atom = !strcmp("H", the_coords[i].name);
60			}
61
62	gezelter	864	if(!skip_atom){
63
64	tim	2750	gb_atom = !strcmp("GB", the_coords[i].name);
65			gbdp_atom = !strcmp("GBDP", the_coords[i].name);
66
67			if (gb_atom) {
68			fprintf(out_file,
69			"make_%s_ellipse( %lf, %lf, %lf, %lf, %lf, %lf, %lf)\n",
70			the_coords[i].name,
71			the_coords[i].x,
72			the_coords[i].z,
73			the_coords[i].y,
74			the_coords[i].charge,
75			the_coords[i].ux,
76			the_coords[i].uz,
77			the_coords[i].uy);
78			} else {
79			if (gbdp_atom) {
80			fprintf(out_file,
81			"make_%s_shaded_ellipse( %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf)\n",
82			the_coords[i].name,
83			the_coords[i].x,
84			the_coords[i].z,
85			the_coords[i].y,
86			the_coords[i].charge,
87			the_coords[i].ux,
88			the_coords[i].uz,
89			the_coords[i].uy,
90			the_coords[i].vx,
91			the_coords[i].vz,
92			the_coords[i].vy);
93			} else {
94			fprintf(out_file,
95			"make_%s_atom( %lf, %lf, %lf )\n",
96			the_coords[i].name,
97			the_coords[i].x,
98			the_coords[i].z,
99			the_coords[i].y);
100			}
101			}
102	mmeineke	60	}
103			}
104	tim	2750
105	mmeineke	60	fprintf(out_file,
106			"\n"
107			"\n");
108	gezelter	864	}
109
110
111			if (d_vectors) {
112
113			fprintf(out_file,
114			"//************************************************************\n"
115			"// The list of vectors\n"
116			"//************************************************************\n"
117			"\n"
118			"\n");
119	mmeineke	60
120	gezelter	864	for(i = 0; i < n_atoms; i++){
121
122			if (the_coords[i].hasVector) {
123			fprintf(out_file,
124			"make_%s_vector(%lf, %lf, %lf, %lf, %lf, %lf)\n",
125			the_coords[i].name,
126			the_coords[i].x,
127			the_coords[i].z,
128			the_coords[i].y,
129			the_coords[i].ux,
130			the_coords[i].uz,
131			the_coords[i].uy);
132			}
133			}
134
135			fprintf(out_file,
136			"\n"
137			"\n");
138	mmeineke	60	}
139	gezelter	864
140	mmeineke	60	if(d_bonds){
141
142			fprintf(out_file,
143			"//************************************************************\n"
144			"// The list of bonds\n"
145			"//************************************************************\n"
146			"\n"
147			"\n");
148
149	mmeineke	508	if( bl_head == NULL ) make_bonds(the_coords, n_atoms);
150	mmeineke	60
151			current_bond = bl_head->next;
152
153			while(current_bond != NULL){
154
155			skip_bond = 0;
156
157			i = current_bond->the_bond.i;
158			j = current_bond->the_bond.j;
159
160			if(!d_hydrogens){
161
162			test1 = !strcmp("H", the_coords[i].name);
163			test2 = !strcmp("H", the_coords[j].name);
164
165			skip_bond = (test1 \|\| test2);
166			}
167
168			if(!skip_bond){
169
170			dx = (the_coords[j].x - the_coords[i].x) / 2.0;
171			dy = (the_coords[j].y - the_coords[i].y) / 2.0;
172			dz = (the_coords[j].z - the_coords[i].z) / 2.0;
173
174			fprintf(out_file,
175			"make_%s_bond( %lf, %lf, %lf, %lf, %lf, %lf )\n",
176			the_coords[i].name,
177			the_coords[i].x,
178			the_coords[i].z,
179			the_coords[i].y,
180			(the_coords[i].x + dx),
181			(the_coords[i].z + dz),
182			(the_coords[i].y + dy));
183
184			fprintf(out_file,
185			"make_%s_bond( %lf, %lf, %lf, %lf, %lf, %lf )\n",
186			the_coords[j].name,
187			the_coords[j].x,
188			the_coords[j].z,
189			the_coords[j].y,
190			(the_coords[j].x - dx),
191			(the_coords[j].z - dz),
192			(the_coords[j].y - dy));
193
194			fprintf(out_file, "\n");
195			}
196
197			current_bond = current_bond->next;
198			}
199
200	mmeineke	508	if( regenerateBonds )clean_bonds();
201	mmeineke	60	}
202			}
203
204
205			void make_bonds(struct coords *the_coords, int n_atoms){
206
207			int i, j; /counters /
208			struct linked_bond_list temp_bond; /bond place holder */
209
210			const double bond_fudge = 1.12; // a fudge factor
211			struct atom type_i, type_j; /* holds the atom types */
212
213			int test; /* booleans */
214			double dx, dy, dz, dr2, dcv, dcv2; // used to determine bond existence
215
216
217			bl_head = (struct linked_bond_list *)malloc(sizeof(struct linked_bond_list));
218			bl_head->next = NULL;
219
220			for(i = 0; i < (n_atoms - 1); i++){
221
222			for(j = (i+1); j < n_atoms; j++){
223
224			dx = the_coords[j].x - the_coords[i].x;
225			dy = the_coords[j].y - the_coords[i].y;
226			dz = the_coords[j].z - the_coords[i].z;
227
228			dr2 = dx * dx + dy * dy + dz * dz;
229
230			test = !findAtomType(the_coords[i].name, &type_i);
231			if(test){
232			fprintf(stderr, "Atom Type %s, not found!\n",
233			the_coords[i].name);
234			exit(8);
235			}
236
237			test = !findAtomType(the_coords[j].name, &type_j);
238			if(test){
239			fprintf(stderr, "Atom Type %s, not found!\n",
240			the_coords[j].name);
241			exit(8);
242			}
243
244
245			dcv = bond_fudge * (type_i.covalentRadii + type_j.covalentRadii);
246			dcv2 = dcv * dcv;
247
248			if(dr2 <= dcv2){
249
250			temp_bond =
251			(struct linked_bond_list *)malloc(sizeof(struct linked_bond_list));
252
253			bl_head->the_bond.i = i;
254			bl_head->the_bond.j = j;
255
256			temp_bond->next = bl_head;
257			bl_head = temp_bond;
258			}
259			}
260			}
261			}
262
263
264			void clean_bonds(){
265			struct linked_bond_list *current_bond;
266			struct linked_bond_list next_bond; / place holders */
267
268
269			current_bond = bl_head->next;
270
271			while(current_bond != NULL){
272
273			next_bond = current_bond->next;
274			free(current_bond);
275			current_bond = next_bond;
276			}
277
278			bl_head->next = NULL;
279	mmeineke	508	free( bl_head );
280			bl_head = NULL;
281	mmeineke	60	}
282
283
284			void make_header_macros(FILE *out_file){
285
286			struct linked_atom *type_list; // list of all atom types
287			struct linked_atom *current_type; // current atom type
288
289			char *name;
290			double red, green, blue;
291			double radius;
292
293			type_list = get_type_list();
294			current_type = type_list->next;
295
296			while(current_type != NULL){
297
298			name = current_type->myAtom.name;
299			radius = current_type->myAtom.vanDerWallRadii;
300			red = ((double)current_type->myAtom.red) / 255.0;
301			green = ((double)current_type->myAtom.green) / 255.0;
302			blue = ((double)current_type->myAtom.blue) / 255.0;
303
304
305
306			fprintf(out_file,
307			"//****************************************************\n"
308			"// DEFINE %s MACROS\n"
309			"//****************************************************\n"
310			"\n"
311			"#macro make_%s_bond "
312			"(end_1x, end_1y, end_1z, end_2x, end_2y, end_2z)\n"
313			"\n"
314			" #local x1 = end_1x;\n"
315			" #local y1 = end_1y;\n"
316			" #local z1 = end_1z;\n"
317			" #local x2 = end_2x;\n"
318			" #local y2 = end_2y;\n"
319			" #local z2 = end_2z;\n"
320			"\n"
321			" #if(ROTATE)\n"
322	gezelter	868	" #local x1_new = rotatePointX + A11 * (x1-rotatePointX) + A12 * (y1-rotatePointY) + A13 * (z1-rotatePointZ);\n"
323			" #local y1_new = rotatePointY + A21 * (x1-rotatePointX) + A22 * (y1-rotatePointY) + A23 * (z1-rotatePointZ);\n"
324			" #local z1_new = rotatePointZ + A31 * (x1-rotatePointX) + A32 * (y1-rotatePointY) + A33 * (z1-rotatePointZ);\n"
325	mmeineke	60	"\n"
326	gezelter	868	" #local x2_new = rotatePointX + A11 * (x2-rotatePointX) + A12 * (y2-rotatePointY) + A13 * (z2-rotatePointZ);\n"
327			" #local y2_new = rotatePointY + A21 * (x2-rotatePointX) + A22 * (y2-rotatePointY) + A23 * (z2-rotatePointZ);\n"
328			" #local z2_new = rotatePointZ + A31 * (x2-rotatePointX) + A32 * (y2-rotatePointY) + A33 * (z2-rotatePointZ);\n"
329	mmeineke	60	"\n"
330			" #else\n"
331			" #local x1_new = x1;"
332			" #local y1_new = y1;"
333			" #local z1_new = z1;"
334			"\n"
335			" #local x2_new = x2;"
336			" #local y2_new = y2;"
337			" #local z2_new = z2;"
338			"\n"
339			" #end\n"
340			"\n"
341			" cylinder{\n"
342			" < x1_new, y1_new, z1_new >,\n"
343			" < x2_new, y2_new, z2_new >,\n"
344			" BOND_RADIUS\n"
345			" texture{\n"
346			" pigment{ rgb < %lf, %lf, %lf > }\n"
347			" finish{\n"
348			" ambient .2\n"
349			" diffuse .6\n"
350			" specular 1\n"
351			" roughness .001\n"
352			" metallic\n"
353			" }\n"
354			" }\n"
355			" }\n"
356			"#end\n"
357			"#macro make_%s_atom "
358			"(center_x, center_y, center_z)\n"
359			"\n"
360			" #local x1 = center_x;\n"
361			" #local y1 = center_y;\n"
362			" #local z1 = center_z;\n"
363			"\n"
364			" #if(ROTATE)\n"
365			"\n"
366	gezelter	868	" #local x1_new = rotatePointX + A11 * (x1-rotatePointX) + A12 * (y1-rotatePointY) + A13 * (z1-rotatePointZ);\n"
367			" #local y1_new = rotatePointY + A21 * (x1-rotatePointX) + A22 * (y1-rotatePointY) + A23 * (z1-rotatePointZ);\n"
368			" #local z1_new = rotatePointZ + A31 * (x1-rotatePointX) + A32 * (y1-rotatePointY) + A33 * (z1-rotatePointZ);\n"
369	mmeineke	60	"\n"
370			" #else\n"
371			"\n"
372			" #local x1_new = x1;"
373			" #local y1_new = y1;"
374			" #local z1_new = z1;"
375			"\n"
376			" #end\n"
377			"\n"
378			" sphere{\n"
379			" < x1_new, y1_new, z1_new >,\n"
380			" ATOM_SPHERE_FACTOR * %lf\n"
381			" texture{\n"
382			" pigment{ rgb < %lf, %lf, %lf > }\n"
383			" finish{\n"
384			" ambient .2\n"
385			" diffuse .6\n"
386			" specular 1\n"
387			" roughness .001\n"
388			" metallic\n"
389			" }\n"
390			" }\n"
391			" }\n"
392			"#end\n"
393	tim	2750	"#macro make_%s_ellipse "
394			"(center_x, center_y, center_z, ecc, u_x, u_y, u_z)\n"
395			"\n"
396			" #local x1 = center_x;\n"
397			" #local y1 = center_y;\n"
398			" #local z1 = center_z;\n"
399			" #local x2 = u_x;\n"
400			" #local y2 = u_y;\n"
401			" #local z2 = u_z;\n"
402			"\n"
403			" #if(ROTATE)\n"
404			"\n"
405			" #local x1_new = rotatePointX + A11 * (x1-rotatePointX) + A12 * (y1-rotatePointY) + A13 * (z1-rotatePointZ);\n"
406			" #local y1_new = rotatePointY + A21 * (x1-rotatePointX) + A22 * (y1-rotatePointY) + A23 * (z1-rotatePointZ);\n"
407			" #local z1_new = rotatePointZ + A31 * (x1-rotatePointX) + A32 * (y1-rotatePointY) + A33 * (z1-rotatePointZ);\n"
408			"\n"
409			" #local x2_new = rotatePointX + A11 * (x2-rotatePointX) + A12 * (y2-rotatePointY) + A13 * (z2-rotatePointZ);\n"
410			" #local y2_new = rotatePointY + A21 * (x2-rotatePointX) + A22 * (y2-rotatePointY) + A23 * (z2-rotatePointZ);\n"
411			" #local z2_new = rotatePointZ + A31 * (x2-rotatePointX) + A32 * (y2-rotatePointY) + A33 * (z2-rotatePointZ);\n"
412			"\n"
413			" #else\n"
414			"\n"
415			" #local x1_new = x1;"
416			" #local y1_new = y1;"
417			" #local z1_new = z1;"
418			"\n"
419			" #local x2_new = x2;"
420			" #local y2_new = y2;"
421			" #local z2_new = z2;"
422			"\n"
423			" #end\n"
424			"\n"
425			" #local myUlen = sqrt(x2_newx2_new + y2_newy2_new + z2_new*z2_new);\n"
426			" #local uux = x2_new / myUlen;\n"
427			" #local uuy = y2_new / myUlen;\n"
428			" #local uuz = z2_new / myUlen;\n"
429			" #local myTheta = -degrees(acos(uuz));\n"
430			" #local myPsi = -degrees(atan(uux/uuy));\n"
431			" #local myScale = ATOM_SPHERE_FACTOR * %lf;\n"
432			"\n"
433			" sphere{\n"
434			" < 0, 0, 0 >, 1\n"
435			" texture{\n"
436			" pigment{\n"
437			" average\n"
438			" pigment_map{\n"
439			" [1.0 grad1]\n"
440			" [1.0 grad2]\n"
441			" [1.0 grad3]\n"
442			" [5.0 gradz]\n"
443			" }\n"
444			" }\n"
445			" finish{\n"
446			" ambient .2\n"
447			" diffuse .6\n"
448			" specular 1\n"
449			" roughness .001\n"
450			" metallic\n"
451			" }\n"
452			" }\n"
453			" scale<myScale,myScale,ecc*myScale>\n"
454			" rotate<myTheta,0,myPsi>\n"
455			" translate< x1_new, y1_new, z1_new>\n"
456			" }\n"
457			"#end\n"
458			"#macro make_%s_shaded_ellipse "
459			"(center_x, center_y, center_z, ecc, u_x, u_y, u_z, v_x, v_y, v_z)\n"
460			"\n"
461			" #local x1 = center_x;\n"
462			" #local y1 = center_y;\n"
463			" #local z1 = center_z;\n"
464			" #local x2 = u_x;\n"
465			" #local y2 = u_y;\n"
466			" #local z2 = u_z;\n"
467			" #local x3 = v_x;\n"
468			" #local y3 = v_y;\n"
469			" #local z3 = v_z;\n"
470			"\n"
471			" #if(ROTATE)\n"
472			"\n"
473			" #local x1_new = rotatePointX + A11 * (x1-rotatePointX) + A12 * (y1-rotatePointY) + A13 * (z1-rotatePointZ);\n"
474			" #local y1_new = rotatePointY + A21 * (x1-rotatePointX) + A22 * (y1-rotatePointY) + A23 * (z1-rotatePointZ);\n"
475			" #local z1_new = rotatePointZ + A31 * (x1-rotatePointX) + A32 * (y1-rotatePointY) + A33 * (z1-rotatePointZ);\n"
476			"\n"
477			" #local x2_new = rotatePointX + A11 * (x2-rotatePointX) + A12 * (y2-rotatePointY) + A13 * (z2-rotatePointZ);\n"
478			" #local y2_new = rotatePointY + A21 * (x2-rotatePointX) + A22 * (y2-rotatePointY) + A23 * (z2-rotatePointZ);\n"
479			" #local z2_new = rotatePointZ + A31 * (x2-rotatePointX) + A32 * (y2-rotatePointY) + A33 * (z2-rotatePointZ);\n"
480			"\n"
481			" #local x3_new = rotatePointX + A11 * (x3-rotatePointX) + A12 * (y3-rotatePointY) + A13 * (z3-rotatePointZ);\n"
482			" #local y3_new = rotatePointY + A21 * (x3-rotatePointX) + A22 * (y3-rotatePointY) + A23 * (z3-rotatePointZ);\n"
483			" #local z3_new = rotatePointZ + A31 * (x3-rotatePointX) + A32 * (y3-rotatePointY) + A33 * (z3-rotatePointZ);\n"
484			"\n"
485			" #else\n"
486			"\n"
487			" #local x1_new = x1;"
488			" #local y1_new = y1;"
489			" #local z1_new = z1;"
490			"\n"
491			" #local x2_new = x2;"
492			" #local y2_new = y2;"
493			" #local z2_new = z2;"
494			"\n"
495			" #local x3_new = x3;"
496			" #local y3_new = y3;"
497			" #local z3_new = z3;"
498			"\n"
499			" #end\n"
500			"\n"
501			" #local myUlen = sqrt(x2_newx2_new + y2_newy2_new + z2_new*z2_new);\n"
502			" #local uux = x2_new / myUlen;\n"
503			" #local uuy = y2_new / myUlen;\n"
504			" #local uuz = z2_new / myUlen;\n"
505			" #local myVlen = sqrt(x3_newx3_new + y3_newy3_new + z3_new*z3_new);\n"
506			" #local vvx = x3_new / myVlen;\n"
507			" #local vvy = y3_new / myVlen;\n"
508			" #local vvz = z3_new / myVlen;\n"
509			"\n"
510			" #local myTheta = degrees(acos(uuz));\n"
511			" #local myPsi = -degrees(atan(uux/uuy));\n"
512			" #local myPhi = degrees(acos(vvz));\n"
513			" #local myScale = ATOM_SPHERE_FACTOR * %lf;\n"
514			"\n"
515			" sphere{\n"
516			" < 0, 0, 0 >, 1\n"
517			" texture{\n"
518			" pigment{\n"
519			" average\n"
520			" pigment_map{\n"
521			" [1.0 grad1]\n"
522			" [1.0 grad2]\n"
523			" [1.0 grad3]\n"
524			" [5.0 gradz]\n"
525			" }\n"
526			" }\n"
527			" finish{\n"
528			" ambient .2\n"
529			" diffuse .6\n"
530			" specular 1\n"
531			" roughness .001\n"
532			" metallic\n"
533			" }\n"
534			" }\n"
535			" scale<myScale,myScale,ecc*myScale>\n"
536			" rotate<myTheta,myPhi,myPsi>\n"
537			" translate< x1_new, y1_new, z1_new>\n"
538			" }\n"
539			"#end\n"
540	gezelter	864	"#macro make_%s_vector "
541			"(center_x, center_y, center_z, ux, uy, uz)\n"
542	mmeineke	60	"\n"
543	gezelter	864	" #local vx = VECTOR_SCALE * ux;\n"
544			" #local vy = VECTOR_SCALE * uy;\n"
545			" #local vz = VECTOR_SCALE * uz;\n"
546			" #local x1 = center_x - 0.5 * vx;\n"
547			" #local y1 = center_y - 0.5 * vy;\n"
548			" #local z1 = center_z - 0.5 * vz;\n"
549			" #local x2 = center_x + 0.5 * vx;\n"
550			" #local y2 = center_y + 0.5 * vy;\n"
551			" #local z2 = center_z + 0.5 * vz;\n"
552			" #local v2 = vxvx + vyvy + vz*vz;\n"
553			" #local vl = sqrt(v2);\n"
554			" #local x3 = x1 + vx * (1.0 - CONE_FRACTION);\n"
555			" #local y3 = y1 + vy * (1.0 - CONE_FRACTION);\n"
556			" #local z3 = z1 + vz * (1.0 - CONE_FRACTION);\n"
557			"\n"
558			" #if(ROTATE)\n"
559	gezelter	868	" #local x1_new = rotatePointX + A11 * (x1-rotatePointX) + A12 * (y1-rotatePointY) + A13 * (z1-rotatePointZ);\n"
560			" #local y1_new = rotatePointY + A21 * (x1-rotatePointX) + A22 * (y1-rotatePointY) + A23 * (z1-rotatePointZ);\n"
561			" #local z1_new = rotatePointZ + A31 * (x1-rotatePointX) + A32 * (y1-rotatePointY) + A33 * (z1-rotatePointZ);\n"
562	gezelter	864	"\n"
563	gezelter	868	" #local x2_new = rotatePointX + A11 * (x2-rotatePointX) + A12 * (y2-rotatePointY) + A13 * (z2-rotatePointZ);\n"
564			" #local y2_new = rotatePointY + A21 * (x2-rotatePointX) + A22 * (y2-rotatePointY) + A23 * (z2-rotatePointZ);\n"
565			" #local z2_new = rotatePointZ + A31 * (x2-rotatePointX) + A32 * (y2-rotatePointY) + A33 * (z2-rotatePointZ);\n"
566	gezelter	864	"\n"
567	gezelter	868	" #local x3_new = rotatePointX + A11 * (x3-rotatePointX) + A12 * (y3-rotatePointY) + A13 * (z3-rotatePointZ);\n"
568			" #local y3_new = rotatePointY + A21 * (x3-rotatePointX) + A22 * (y3-rotatePointY) + A23 * (z3-rotatePointZ);\n"
569			" #local z3_new = rotatePointZ + A31 * (x3-rotatePointX) + A32 * (y3-rotatePointY) + A33 * (z3-rotatePointZ);\n"
570	gezelter	864	"\n"
571			" #else\n"
572			" #local x1_new = x1;"
573			" #local y1_new = y1;"
574			" #local z1_new = z1;"
575			"\n"
576			" #local x2_new = x2;"
577			" #local y2_new = y2;"
578			" #local z2_new = z2;"
579			"\n"
580			" #local x3_new = x3;"
581			" #local y3_new = y3;"
582			" #local z3_new = z3;"
583			"\n"
584			" #end\n"
585			"\n"
586			" cylinder{\n"
587			" < x1_new, y1_new, z1_new >,\n"
588			" < x3_new, y3_new, z3_new >,\n"
589			" STICK_RADIUS\n"
590			" texture{\n"
591			" pigment{ rgb < %lf, %lf, %lf > }\n"
592			" finish{\n"
593			" ambient .2\n"
594			" diffuse .6\n"
595			" specular 1\n"
596			" roughness .001\n"
597			" metallic\n"
598			" }\n"
599			" }\n"
600			" }\n"
601			" cone{\n"
602			" < x2_new, y2_new, z2_new >, 0.0\n"
603			" < x3_new, y3_new, z3_new >, CONE_RADIUS\n"
604			" texture{\n"
605			" pigment{ rgb < %lf, %lf, %lf > }\n"
606			" finish{\n"
607			" ambient .2\n"
608			" diffuse .6\n"
609			" specular 1\n"
610			" roughness .001\n"
611			" metallic\n"
612			" }\n"
613			" }\n"
614			" }\n"
615			"#end\n"
616			"\n"
617	mmeineke	60	"\n",
618			name,
619			name,
620			red, green, blue,
621			name,
622			radius,
623	gezelter	864	red, green, blue,
624			name,
625	tim	2750	radius,
626			name,
627			radius,
628			name,
629	gezelter	864	red, green, blue,
630	mmeineke	60	red, green, blue);
631
632			current_type = current_type->next;
633			}
634			}