cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole_species.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*CO_Init called from cdInit to initialize co routines */
4 /*CO_update_chem_rates update rate coefficients, only temp part - in mole_co_etc.c */
5 #include "cdstd.h"
6 #include "cddefines.h"
7 #include "deuterium.h"
8 #include "elementnames.h"
9 #include "h2.h"
10 #include "hmi.h"
11 #include "dense.h"
12 #include "grainvar.h"
13 #include "trace.h"
14 #include "abund.h"
15 #include "mole_priv.h"
16 #include "mole.h"
17 #include "parse_species.h"
18 /*lint -e778 constant expression evaluates to 0 in operation '-' */
19 
20 /* CO_update_chem_rates update rate coefficients, only temp part - in
21  * mole_co_etc.c called in conv_base before any chemistry or
22  * ionization is done */
23 
25 
26 STATIC void read_species_file( string filename, bool lgCreateIsotopologues = true );
27 STATIC void newelement(const char label[], int ipion);
28 STATIC void newisotope( const count_ptr<chem_element> &el, int massNumberA,
29  realnum mass_amu, double frac );
31 // newspecies is overloaded. The first one just calls the second with the additional lgCreateIsotopologues set true
32 STATIC molecule *newspecies(const char label[], enum spectype type,
33  enum mole_state state, realnum form_enthalpy);
34 STATIC molecule *newspecies(const char label[], enum spectype type,
35  enum mole_state state, realnum form_enthalpy, bool lgCreateIsotopologues);
36 STATIC bool isactive(const molecule &mol);
37 STATIC bool ispassive(const molecule &mol);
38 STATIC void SetIsotopeFractions( const vector<bool>& lgResolveNelem );
39 
40 namespace mole_priv {
41  map <string,count_ptr<molecule> > spectab;
42  map <string,count_ptr<mole_reaction> > reactab;
43  map <string,count_ptr<chem_element> > elemtab;
44  map <string,count_ptr<mole_reaction> > functab;
45 }
46 
47 map <string,size_t > nuclidetab;
48 typedef map<string,size_t >::iterator chem_nuclide_i;
49 
56 vector<molecule *> groupspecies;
57 
58 
59 namespace
60 {
61  class MoleCmp : public binary_function<const count_ptr<molecule>,
62  const count_ptr<molecule>,bool>
63  {
64  public:
65  bool operator()(const count_ptr<molecule> &mol1,
66  const count_ptr<molecule> &mol2) const
67  {
68  return mol1->compare(*mol2) < 0;
69  }
70  bool operator()(const molecule *mol1, const molecule *mol2) const
71  {
72  return mol1->compare(*mol2) < 0;
73  }
74  };
75 }
76 
77 void t_mole_global::sort(t_mole_global::MoleculeList::iterator start,
78  t_mole_global::MoleculeList::iterator end)
79 {
80  std::sort(start,end,MoleCmp());
81 }
83 {
84  std::sort(start,end,MoleCmp());
85 }
86 
87 STATIC void SetIsotopeFractions( const vector<bool>& lgResolveNelem )
88 {
89  DEBUG_ENTRY( "SetIsotopeFractions()" );
90 
91  for( int ielm = 0; ielm < LIMELM; ielm++ )
92  {
93  long Z = ielm;
94  for( int j = 0; j < abund.IsoAbn[ielm].getNiso(); j++ )
95  {
96  long A = abund.IsoAbn[ielm].getAiso( j );
97  realnum mass = abund.IsoAbn[ielm].getMass( A );
98  realnum frac = abund.IsoAbn[ielm].getAbn( A );
99 
100  if( (unsigned)Z <= lgResolveNelem.size() && lgResolveNelem[Z] )
101  newisotope( element_list[Z], A, mass, frac );
102  // always do this to continue history of predicting 13CO
103  else if( Z==ipCARBON )
104  newisotope( element_list[Z], A, mass, frac );
105  }
106  }
107 
108  return;
109 }
110 
112 {
113  DEBUG_ENTRY( "mole_global::make_species()" );
114 
115  long int i;
116  molecule *sp;
117 
118  null_element = new chem_element(0,"") ;
119  null_nuclide = new( chem_nuclide );
120  null_molezone = new( molezone );
121  null_molezone->den = 0.;
122 
123  /* set up concordance of elemental species to external Cloudy indices */
124  for (i=0;i<LIMELM;i++)
125  {
127  }
128 
129  // flip this to treat deuterium
130  if( deut.lgElmtOn )
131  lgTreatIsotopes[ipHYDROGEN] = true;
132 
133  // read and define isotopes, set default fractionations
135 
136 
138  {
139  SetDeuteriumFractionation( element_list[ipHYDROGEN]->isotopes[2]->frac );
140  SetGasPhaseDeuterium( dense.gas_phase[ipHYDROGEN] );
142  }
143 
144  for( long nelem=0; nelem<LIMELM; ++nelem )
145  {
146  realnum mass_amu = MeanMassOfElement( element_list[nelem] );
147  //define generic mean-abundance isotopes
148  newisotope( element_list[nelem], -1, mass_amu, 1.0 );
149  }
150 
151  /* set up properties of molecular species -- chemical formulae,
152  array indices, elementary components (parsed from formula),
153  status within CO network, location of stored value external
154  to CO network (must be floating point). */
155 
156  /* Sizes of different parts of network are calculated by increments in newspecies */
158  /* Enthalpies of formation taken from
159  * >>refer mol Le Teuff, Y. E., Millar, T. J., & Markwick, A. J.,2000, A&AS, 146, 157
160  */
161 
162  /* Zero density pseudo-species to return when molecule is switched off */
163  null_mole = newspecies(" ",OTHER,MOLE_NULL, 0.);
164  null_mole->index = -1;
165 
166  read_species_file( "chem_species.dat" );
167 
169  {
170  /* What are formation enthalpies of COgrn, H2Ogrn, OHgrn? For
171  present, take grn as standard state, and neglect adsorption enthalpy.
172 
173  -- check e.g. http://www.arxiv.org/abs/astro-ph/0702322 for CO adsorption energy.
174  */
175  if (0)
176  {
177  read_species_file( "chem_species_grn.dat" );
178  }
179  else
180  {
181  (void) newspecies("COgrn ",MOLECULE,MOLE_ACTIVE, -113.8f);
182  (void) newspecies("H2Ogrn",MOLECULE,MOLE_ACTIVE, -238.9f);
183  (void) newspecies("OHgrn ",MOLECULE,MOLE_ACTIVE, 38.4f);
184  //sp = newspecies("Hgrn ",MOLECULE,MOLE_ACTIVE, 216.f);
185  }
186  }
187 
188  /* Add passive species to complete network */
189  sp = newspecies("e- ",OTHER,MOLE_PASSIVE, 0.0f);
190  sp->charge = -1; sp->mole_mass = (realnum)ELECTRON_MASS; /* Augment properties for this non-molecular species */
191  (void) newspecies("grn ",OTHER,MOLE_PASSIVE, 0.0f);
192  (void) newspecies("PHOTON",OTHER,MOLE_PASSIVE, 0.0f);
193  (void) newspecies("CRPHOT",OTHER,MOLE_PASSIVE, 0.0f);
194  (void) newspecies("CRP ",OTHER,MOLE_PASSIVE, 0.0f);
195 
197  (void) newspecies("H- ",MOLECULE,MOLE_ACTIVE, 143.2f);
199  {
200  (void) newspecies("H2* ",MOLECULE,MOLE_ACTIVE,
201  h2.ENERGY_H2_STAR * KJMOL1CM);
202  }
203 
204  if( deut.lgElmtOn )
205  {
206  read_species_file( "chem_species_deuterium.dat", false );
207  }
208 
209  /* Add species for all other elements and their first ions -- couple to network at least via H- */
210  for( ChemNuclideList::iterator atom = nuclide_list.begin();
211  atom != nuclide_list.end(); ++atom)
212  {
213  if ( !(*atom)->lgHasLinkedIon() )
214  continue;
215  unsigned long int nelem = (*atom)->el->Z-1;
216 
217  for( unsigned long ion=0; ion<=nelem+1; ion++ )
218  {
220  char temp[CHARS_ION_STAGE+1];
221  if( ion==0 )
222  temp[0] = '\0';
223  else if( ion==1 )
224  sprintf( temp, "+" );
225  else if( ion < 100 )
226  sprintf( temp, "+%ld", ion );
227  else
228  TotalInsanity();
229  sprintf( str, "%s%s", (*atom)->label().c_str(), temp );
230  if (findspecies(str) == null_mole)
231  {
232  sp = newspecies(str,MOLECULE,MOLE_ACTIVE, 0.f);
233  fixit("populate these in a local update");
234 #if 0
235  if( sp != NULL )
236  {
237  sp->levels = NULL;
238  sp->numLevels = 0;
239  }
240 #endif
241  }
242  }
243  }
244 
245  return;
246 }
247 
249 {
250  DEBUG_ENTRY( "mole_make_list()" );
251 
252  /* Create linear list of species and populate it... */
254 
255  /* ...first active species */
256  long int i = 0;
257  for (molecule_i p = mole_priv::spectab.begin(); p != mole_priv::spectab.end(); ++p)
258  {
259  if (isactive(*(p->second)))
260  mole_global.list[i++] = p->second;
261  }
262  ASSERT (i == mole_global.num_calc);
263 
264  /* Sort list into a standard ordering */
266 
267  for (molecule_i p = mole_priv::spectab.begin(); p != mole_priv::spectab.end(); ++p)
268  {
269  if (ispassive(*(p->second)))
270  mole_global.list[i++] = p->second;
271  }
272  ASSERT (i == mole_global.num_total);
273 
274  /* Set molecule indices to order of list just created */
275  for(i=0;i<mole_global.num_total;i++)
276  {
277  mole_global.list[i]->index = i;
278  }
279 
280  for(i=0;i<mole_global.num_total;i++)
281  {
282  if( !mole_global.list[i]->parentLabel.empty() )
283  {
284  long parentIndex = findspecies( mole_global.list[i]->parentLabel.c_str() )->index;
285  mole_global.list[i]->parentIndex = parentIndex;
286  }
287  else
288  mole_global.list[i]->parentIndex = -1;
289  }
290 
291  /* Register the atomic ladders */
292  for(i=0;i<mole_global.num_total;i++)
293  {
294  molecule *sp = &(*mole_global.list[i]);
295  if (sp->isMonatomic())
296  {
297  ASSERT( (int)sp->nNuclide.size() == 1 );
298  count_ptr<chem_nuclide> atom = sp->nNuclide.begin()->first;
299  ASSERT(sp->charge <= atom->el->Z);
300  if(sp->charge >= 0 && sp->lgGas_Phase)
301  {
302  atom->ipMl[sp->charge] = i;
303  }
304  }
305  }
306 
307  return;
308 
309 }
310 
311 STATIC void read_species_file( string filename, bool lgCreateIsotopologues )
312 {
313  DEBUG_ENTRY( "read_species_file()" );
314 
315  fstream ioDATA;
316  open_data( ioDATA, filename.c_str(), mode_r );
317  string line;
318 
319  while( getline( ioDATA,line ) )
320  {
321  if( line.empty() )
322  break;
323  if( line[0] == '#' )
324  continue;
325  istringstream iss( line );
326  string species;
327  double formation_enthalpy;
328  iss >> species;
329  iss >> formation_enthalpy;
330  ASSERT( iss.eof() );
331  newspecies( species.c_str(), MOLECULE,MOLE_ACTIVE, formation_enthalpy, lgCreateIsotopologues );
332  //fprintf( ioQQQ, "DEBUGGG read_species_file %s\t%f\n", species.c_str(), formation_enthalpy );
333  }
334 
335  return;
336 }
337 /*lint +e778 constant expression evaluates to 0 in operation '-' */
338 
341  vector< int >& numNuclides,
342  string atom_old,
343  string atom_new,
344  string embellishments,
345  vector<string>& newLabels )
346 {
347  DEBUG_ENTRY( "create_isotopologues()" );
348 
349  fixit("make sure atom_new and atom_old are isotopes");
350  fixit("make sure atom_new is not already present");
351 
352  //for( ChemNuclideList::iterator it = atoms.begin(); it != atoms.end(); ++it )
353  for( unsigned position = 0; position < atoms.size(); ++position )
354  {
355  string newLabel;
356  create_isotopologues_one_position( position, atoms, numNuclides, atom_old, atom_new, embellishments, newLabel );
357  if( !newLabel.empty() )
358  newLabels.push_back( newLabel );
359  }
360 
361  return;
362 }
363 
365  unsigned position,
367  vector< int >& numNuclides,
368  string atom_old,
369  string atom_new,
370  string embellishments,
371  string& newLabel )
372 {
373  DEBUG_ENTRY( "create_isotopologues_one_position()" );
374 
375  stringstream str;
376  if( atoms[position]->label() == atom_old )
377  {
378  for( unsigned i=0; i<position; ++i )
379  {
380  str << atoms[i]->label();
381  if( numNuclides[i]>1 )
382  str << numNuclides[i];
383  }
384 
385  if( numNuclides[position] > 1 )
386  {
387  str << atom_old;
388  if( numNuclides[position] > 2 )
389  str << numNuclides[position]-1;
390  }
391 
392  if( position+1 == atoms.size() )
393  str << atom_new;
394 
395  for( unsigned i=position+1; i<atoms.size(); ++i )
396  {
397  if( i==position+1 )
398  {
399  // remove adjacent duplicates
400  if( atom_new == atoms[i]->label() )
401  {
402  str << atoms[i]->label();
403  ASSERT( numNuclides[i] + 1 > 1 );
404  str << numNuclides[i] + 1;
405  }
406  else
407  {
408  str << atom_new;
409  str << atoms[i]->label();
410  if( numNuclides[i] > 1 )
411  str << numNuclides[i];
412  }
413  }
414  else
415  {
416  str << atoms[i]->label();
417  if( numNuclides[i] > 1 )
418  str << numNuclides[i];
419  }
420  }
421 
422  // add on charge, grn, and excitation embellishments
423  str << embellishments;
424 
425  newLabel = str.str();
426  //fprintf( ioQQQ, "DEBUGGG create_isotopologues_one_position %s\n", newLabel.c_str() );
427  }
428 
429  return;
430 }
431 
432 /* Fill element linking structure */
433 STATIC void newelement(const char label[], int ipion)
434 {
435  char *s;
436 
437  DEBUG_ENTRY( "newelement()" );
438 
439  /* Create private workspace for label; copy and strip trailing whitespace */
440  int len = strlen(label)+1;
441  auto_vec<char> mylab_v(new char[len]);
442  char *mylab = mylab_v.data();
443  strncpy(mylab,label,len);
444  s = strchr(mylab,' ');
445  if (s)
446  *s = '\0';
447 
448  int exists = (mole_priv::elemtab.find(mylab) != mole_priv::elemtab.end());
449  if (!exists)
450  {
451  count_ptr<chem_element> element(new chem_element(ipion+1,mylab));
452  mole_priv::elemtab[element->label] = element;
453  element_list.push_back(element);
454  }
455  return;
456 }
457 
458 /* Fill isotope lists */
459 STATIC void newisotope( const count_ptr<chem_element>& el, int massNumberA,
460  realnum mass_amu, double frac )
461 {
462 
463  DEBUG_ENTRY( "newisotope()" );
464 
465  ASSERT( massNumberA < 3 * LIMELM && ( massNumberA > 0 || massNumberA == -1 ) );
466  ASSERT( mass_amu < 3. * LIMELM && mass_amu > 0. );
467  ASSERT( frac <= 1. + FLT_EPSILON && frac >= 0. );
468 
470  isotope->A = massNumberA;
471  isotope->mass_amu = mass_amu;
472  isotope->frac = frac;
473  isotope->ipMl.resize(el->Z+1);
474  isotope->el = el.get_ptr();
475  for (long int ion = 0; ion < el->Z+1; ion++)
476  isotope->ipMl[ion] = -1; /* Chemical network species indices not yet defined */
477 
478  //int exists = (mole_priv::nuclidetab.find( isotope->label() ) != mole_priv::nuclidetab.end());
479  nuclidetab[ isotope->label() ] = nuclide_list.size();
480  nuclide_list.push_back(isotope);
481  // register with 'parent' element
482  el->isotopes[massNumberA] = isotope;
483 }
484 
486 {
487  DEBUG_ENTRY( "MeanMassOfElement()" );
488 
489  // if no isotopes have been defined, just use the mean mass stored elsewhere
490  if( el->isotopes.size()==0 )
491  return dense.AtomicWeight[el->Z-1];
492 
493  realnum averageMass = 0., fracsum = 0.;
494  for( isotopes_i it = el->isotopes.begin(); it != el->isotopes.end(); ++it )
495  {
496  fracsum += it->second->frac;
497  averageMass += it->second->mass_amu * it->second->frac;
498  }
499  ASSERT( fp_equal( fracsum, realnum(1.) ) );
500 
501  return averageMass;
502 }
503 
505  const char label[],
506  enum spectype type,
507  enum mole_state state,
508  realnum form_enthalpy)
509 {
510  return newspecies( label, type, state, form_enthalpy, true);
511 }
512 
513 /* Parse species string to find constituent atoms, charge etc. */
515  const char label[],
516  enum spectype type,
517  enum mole_state state,
518  realnum form_enthalpy,
519  bool lgCreateIsotopologues )/* formation enthalpy at 0K */
520 {
521  DEBUG_ENTRY( "newspecies()" );
522 
523  int exists;
524  ChemNuclideList nuclidesLeftToRight;
525  vector< int > numNuclides;
526  string embellishments;
527  char *s;
528  count_ptr<molecule> mol(new molecule);
529 
530  mol->parentLabel.clear();
531  mol->isEnabled = true;
532  mol->charge = 0;
533  mol->lgExcit = false;
534  mol->mole_mass = 0.0;
535  mol->state = state;
536  mol->lgGas_Phase = true;
537  mol->form_enthalpy = form_enthalpy;
538  mol->groupnum = -1;
539  mol->n_react = 0;
540 
541  /* Create private workspace for label; copy and strip trailing whitespace */
542  int len = strlen(label)+1;
543  auto_vec<char> mylab_v(new char[len]);
544  char *mylab = mylab_v.data();
545  strncpy(mylab,label,len);
546  s = strchr(mylab,' ');
547  if (s)
548  *s = '\0';
549  mol->label = mylab;
550 
551  if(type == MOLECULE)
552  {
553  if( parse_species_label( mylab, nuclidesLeftToRight, numNuclides, embellishments, mol->lgExcit, mol->charge, mol->lgGas_Phase ) == false )
554  return NULL;
555 
556  for( unsigned i = 0; i < nuclidesLeftToRight.size(); ++i )
557  {
558  mol->nNuclide[ nuclidesLeftToRight[i] ] += numNuclides[i];
559  mol->mole_mass += numNuclides[i] * nuclidesLeftToRight[i]->mass_amu * (realnum)ATOMIC_MASS_UNIT;
560  }
561  }
562 
563  // we also kill H- if molecules are disabled. This is less than ideal,
564  // but physically motivated by the fact that one of the strongest H- sinks
565  // involves formation of H2 (disabled by "no molecules"), while one the
566  // fastest sources is e- recombination (which would still be allowed).
567  if ( (mol->n_nuclei() > 1 || (mol->isMonatomic() && mol->charge==-1) ) && mole_global.lgNoMole)
568  {
569  if( trace.lgTraceMole )
570  fprintf(ioQQQ,"No species %s as molecules off\n",label);
571  return NULL;
572  }
573 
574  if (mol->n_nuclei() > 1 && mole_global.lgNoHeavyMole)
575  {
576  for( nNucs_ri it=mol->nNuclide.rbegin(); it != mol->nNuclide.rend(); --it )
577  {
578  if( it->first->el->Z-1 != ipHYDROGEN )
579  {
580  ASSERT( it->second > 0 );
581  if( trace.lgTraceMole )
582  fprintf(ioQQQ,"No species %s as heavy molecules off\n",label);
583  return NULL;
584  }
585  }
586  }
587 
588  if (speciesOff(mol->label))
589  {
590  fprintf( ioQQQ,"Species %s disabled\n",mol->label.c_str() );
591  return NULL;
592  }
593 
594  /* Insert species into hash table */
595  exists = (mole_priv::spectab.find(mol->label) != mole_priv::spectab.end());
596  if( exists )
597  {
598  fprintf( ioQQQ,"Warning: duplicate species %s - using first one\n",
599  mol->label.c_str() );
600  return NULL;
601  }
602  else
603  mole_priv::spectab[mol->label] = mol;
604 
605  // all map entries should have strictly positive number of nuclei
606  for( nNucs_i it=mol->nNuclide.begin(); it != mol->nNuclide.end(); ++it )
607  ASSERT( it->second > 0 );
608 
609  if (state != MOLE_NULL)
611  if (state == MOLE_ACTIVE)
613 
614  // this is a special case to always treat 13CO (as has long been done)
615  if( lgCreateIsotopologues && type==MOLECULE && mol->label.compare("CO")==0 )
616  {
617  molecule *sp = newspecies( "^13CO", MOLECULE, mol->state, mol->form_enthalpy, false );
618  sp->parentLabel = mol->label;
619  }
620 
621  // create singly-substituted isotopologues
622  if( lgCreateIsotopologues && type==MOLECULE && !mol->isMonatomic() )
623  {
624  for( nNucs_i it1 = mol->nNuclide.begin(); it1 != mol->nNuclide.end(); ++it1 )
625  {
626  for( map<int, count_ptr<chem_nuclide> >::iterator it2 = it1->first->el->isotopes.begin();
627  it2 != it1->first->el->isotopes.end(); ++it2 )
628  {
629  // we don't want to create ^1H isotopologues (only substitute D for H)
630  if( it1->first->el->Z-1 == ipHYDROGEN && it2->second->A != 2 )
631  continue;
632  if( !mole_global.lgTreatIsotopes[it1->first->el->Z-1] )
633  continue;
634  if( it2->second->lgMeanAbundance() )
635  continue;
636  vector<string> isoLabs;
637  create_isotopologues( nuclidesLeftToRight, numNuclides, it1->first->label(), it2->second->label(), embellishments, isoLabs );
638  //fprintf( ioQQQ, " DEBUGGG %10s isotopologues of %10s:", it1->first->label().c_str(), mol->label.c_str() );
639  //for( vector<string>::iterator lab = isoLabs.begin(); lab != isoLabs.end(); ++ lab )
640  // fprintf( ioQQQ, " %10s", lab->c_str() );
641  //fprintf( ioQQQ, "\n" );
642  for( vector<string>::iterator newLabel = isoLabs.begin(); newLabel != isoLabs.end(); ++newLabel )
643  {
644  molecule *sp = newspecies( newLabel->c_str(), MOLECULE, mol->state, mol->form_enthalpy, false );
645  // D is special -- don't set parentLabel
646  if( sp!=NULL && it2->second->A != 2 )
647  sp->parentLabel = mol->label;
648  }
649  }
650  }
651  }
652 
653  return &(*mol);
654 }
655 bool parse_species_label( const char label[], ChemNuclideList &nuclidesLeftToRight, vector<int> &numNuclides, string &embellishments )
656 {
657  bool lgExcit, lgGas_Phase;
658  int charge;
659  bool lgOK = parse_species_label( label, nuclidesLeftToRight, numNuclides, embellishments, lgExcit, charge, lgGas_Phase );
660  return lgOK;
661 }
662 bool parse_species_label( const char label[], ChemNuclideList &nuclidesLeftToRight, vector<int> &numNuclides, string &embellishments,
663  bool &lgExcit, int &charge, bool &lgGas_Phase )
664 {
665  long int i, n, ipAtom;
666  char thisAtom[CHARS_ISOTOPE_SYM];
668  const char *s;
669  int iend = strlen(label);
670 
671  /* Excitation... */
672  s = strpbrk(label,"*");
673  if(s)
674  {
675  lgExcit = true;
676  embellishments = s;
677  iend = s-label;
678  }
679 
680  /* ...Charge */
681  s = strpbrk(label,"+-");
682  if(s)
683  {
684  if(isdigit(*(s+1)))
685  n = atoi(s+1);
686  else
687  n = 1;
688  if(*s == '+')
689  charge = n;
690  else
691  charge = -n;
692  embellishments = s + embellishments;
693  iend = s-label;
694  }
695  /* ...Grain */
696  s = strstr(label,"grn");
697  if(s)
698  {
699  lgGas_Phase = false;
700  embellishments = s + embellishments;
701  iend = s-label;
702  }
703  else
704  {
705  lgGas_Phase = true;
706  }
707  //fprintf( ioQQQ, "DEBUGGG parse_species_label %s\t%d\t%s\t%d\t%s\t\n",
708  // label, (int)strlen(label), label, iend, embellishments.c_str() );
709 
710  /* Now analyse chemical formula */
711  i = 0;
712  while (i < iend && label[i] != ' ' && label[i] != '*')
713  {
714  ipAtom = 0;
715  /* Select next nuclide in species, matches regexp (^\d+)?[A-Z][a-z]? */
716  /* look for isotope prefix */
717  if(label[i]=='^')
718  {
719  thisAtom[ipAtom++] = label[i++];
720  ASSERT( isdigit(label[i]) );
721  thisAtom[ipAtom++] = label[i++];
722  if(isdigit(label[i]))
723  {
724  thisAtom[ipAtom++] = label[i++];
725  }
726  }
727  // should be first character of an element symbol
728  thisAtom[ipAtom++] = label[i++];
729  if( i<iend && islower(label[i]))
730  {
731  thisAtom[ipAtom++] = label[i++];
732  }
733  thisAtom[ipAtom] = '\0';
734  ASSERT(ipAtom <= CHARS_ISOTOPE_SYM);
735 
736  atom = findnuclide(thisAtom);
737  if(atom.get_ptr() == NULL)
738  {
739  fprintf(stderr,"Did not recognize atom at %s in \"%s \"[%ld]\n",thisAtom,label,i);
740  exit(-1);
741  }
742  if(!dense.lgElmtOn[atom->el->Z-1])
743  {
744  if( trace.lgTraceMole )
745  fprintf(ioQQQ,"No species %s as element %s off\n",label,atom->el->label.c_str() );
746  return false;
747  }
748 
749  if(i < iend && isdigit(label[i])) /* If there is >1 of this atom */
750  {
751  n = 0;
752  do {
753  n = 10*n+(long int)(label[i]-'0');
754  i++;
755  } while (i < iend && isdigit(label[i]));
756  }
757  else
758  {
759  n = 1;
760  }
761  nuclidesLeftToRight.push_back( atom );
762  numNuclides.push_back( n );
763  }
764 
765  return true;
766 }
767 STATIC bool isactive(const molecule &mol)
768 {
769  DEBUG_ENTRY( "isactive()" );
770  return mol.state == MOLE_ACTIVE;
771 }
772 STATIC bool ispassive(const molecule &mol)
773 {
774 
775  DEBUG_ENTRY( "ispassive()" );
776  return mol.state == MOLE_PASSIVE;
777 }
778 
779 bool lgDifferByExcitation( const molecule &mol1, const molecule &mol2 )
780 {
781  if( mol1.label == mol2.label + "*" )
782  return true;
783  else if( mol2.label == mol1.label + "*" )
784  return true;
785  else
786  return false;
787 }
788 
789 molecule *findspecies(const char buf[])
790 {
791  DEBUG_ENTRY( "findspecies()" );
792 
793  // strip string of the first space and anything after it
794  string s;
795  for (const char *pb = buf; *pb && *pb != ' '; ++pb)
796  {
797  s += *pb;
798  }
799 
800  const molecule_i p = mole_priv::spectab.find(s);
801 
802  if(p != mole_priv::spectab.end())
803  return &(*p->second);
804  else
805  return null_mole;
806 }
807 
808 molecule *findspecies_validate(const char buf[])
809 {
810  DEBUG_ENTRY( "findspecies_validate()" );
811 
812  // strip string of the first space and anything after it
813  string s;
814  for (const char *pb = buf; *pb && *pb != ' '; ++pb)
815  {
816  s += *pb;
817  }
818 
819  const molecule_i p = mole_priv::spectab.find(s);
820 
821  if(p != mole_priv::spectab.end())
822  return &(*p->second);
823  else
824  {
825  fprintf(ioQQQ," PROBLEM specified species, '%s', is not valid or disabled.\n", s.c_str() );
826  fprintf(ioQQQ," The SAVE SPECIES LABELS command will generate a list of species. Species names are case sensitive.\n");
827  cdEXIT( EXIT_FAILURE );
828  }
829 }
830 
831 molezone *findspecieslocal(const char buf[])
832 {
833  DEBUG_ENTRY( "findspecieslocal()" );
834 
835  molecule *sp = findspecies(buf);
836 
837  if(sp != null_mole)
838  return &mole.species[ sp->index ];
839  else
840  return null_molezone;
841 }
842 
844 {
845  DEBUG_ENTRY( "findspecieslocal_validate()" );
846 
847  molecule *sp = findspecies_validate(buf);
848 
849  return &mole.species[ sp->index ];
850 }
851 
853 {
854  chem_nuclide_i p;
855 
856  DEBUG_ENTRY( "findnuclide()" );
857 
858  p = nuclidetab.find(buf);
859 
860  if(p != nuclidetab.end())
861  return nuclide_list[p->second];
862  else
863  return count_ptr<chem_nuclide>(NULL);
864 }
865 
867 {
868  int i;
869  double den_times_area, den_grains, adsorbed_density;
870 
871  DEBUG_ENTRY( "mole_update_species_cache()" );
872 
873  enum { DEBUG_MOLE = false };
874 
875  /* For rates that are dependent on grain physics. This includes grain density,
876  cross sectional area, and dust temperature of each constituent. Note that
877 
878  gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3
879 
880  is the integrated projected grain surface area per cm^3 of gas for each grain size bin */
881 
882  /* >>chng 06 feb 28, turn off this rate when no grain molecules */
883  /* >>chng 06 dec 05 rjrw: do this in newreact rather than rate */
884  if( gv.lgDustOn() )
885  {
886  den_times_area = den_grains = 0.0;
887  for( size_t nd=0; nd < gv.bin.size(); nd++ )
888  {
889  /* >>chng 06 mar 04, update expression for projected grain surface area, PvH */
890  den_times_area += gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
891  den_grains += gv.bin[nd]->cnv_GR_pCM3;
892  }
893 
894  adsorbed_density = 0.0;
895  for (i=0;i<mole_global.num_total;i++)
896  {
897  if( !mole_global.list[i]->lgGas_Phase && mole_global.list[i]->isIsotopicTotalSpecies() )
898  adsorbed_density += mole.species[i].den;
899  }
900 
901  mole.grain_area = den_times_area;
902  mole.grain_density = den_grains;
903 
904  double mole_cs = 1e-15;
905  if (4*den_times_area <= mole_cs*adsorbed_density)
906  mole.grain_saturation = 1.0;
907  else
908  mole.grain_saturation = ((mole_cs*adsorbed_density)/(4.*den_times_area));
909  }
910  else
911  {
912  mole.grain_area = 0.0;
913  mole.grain_density = 0.0;
914  mole.grain_saturation = 1.0;
915  }
916 
917  if (DEBUG_MOLE)
918  fprintf(ioQQQ,"Dust: %f %f %f\n",
920 
921  for (i=0;i<mole_global.num_total;i++)
922  {
923  if(mole.species[i].location != NULL)
924  {
925  ASSERT( mole_global.list[i]->isIsotopicTotalSpecies() );
926  mole.species[i].den = *(mole.species[i].location);
927  if (DEBUG_MOLE)
928  fprintf(ioQQQ,"Updating %s: %15.8g\n",mole_global.list[i]->label.c_str(),mole.species[i].den);
929  }
930  }
931 
933 }
934 
936 // Finding the total atom density from MoleMap.molElems w
937 {
938  int i;
939 
940  /* These two updates should together maintain the abundance invariant */
941 
942  // Assert invariant
944 
945  // Update total of non-ladder species
948 
949  /* charge on molecules */
950  mole.elec = 0.;
951  for(i=0;i<mole_global.num_calc;i++)
952  {
953  if (mole.species[i].location == NULL && mole_global.list[i]->isIsotopicTotalSpecies())
954  mole.elec += mole.species[i].den*mole_global.list[i]->charge;
955  }
956 
957  // Update ionization ladder species
958  realnum delta = 0.0;
959  long ncpt = 0;
960 
961  for (i=0;i<mole_global.num_total;i++)
962  {
963  if(mole.species[i].location && mole_global.list[i]->state == MOLE_ACTIVE)
964  {
965  realnum new_pop = mole.species[i].den;
966 
967  if( mole_global.list[i]->isMonatomic() )
968  {
969  realnum old_pop = *(mole.species[i].location);
970  long nelem = mole_global.list[i]->nNuclide.begin()->first->el->Z-1;
971  realnum frac = (new_pop-old_pop)/SDIV(new_pop+old_pop+1e-8*
972  dense.gas_phase[nelem]);
973  delta += frac*frac;
974  ++ncpt;
975  }
976 
977  //if( iteration >= 3 && nzone >= 100 )
978  // fprintf( ioQQQ, "DEBUGGG mole_return_ %i\t%s\t%.12e\t%.12e\t%.12e\t%.12e\t%li\n",
979  // i, mole_global.list[i]->label.c_str(), new_pop, old_pop, frac, delta, ncpt );
980  *(mole.species[i].location) = new_pop;
981  }
982  }
983 
984  // Assert invariant
986  return ncpt>0 ? sqrt(delta/ncpt) : 0.f;
987 }
988 
990 {
991  DEBUG_ENTRY( "t_mole_local::set_isotope_abundances()" );
992 
993  // loop over unresolved elements
994  for(ChemNuclideList::iterator atom = nuclide_list.begin(); atom != nuclide_list.end(); ++atom)
995  {
996  if ( !(*atom)->lgMeanAbundance() )
997  continue;
998 
999  // loop over all isotopes of each element
1000  for( isotopes_i it = (*atom)->el->isotopes.begin(); it != (*atom)->el->isotopes.end(); ++it )
1001  {
1002  // skip mean-abundance "isotopes"
1003  if( it->second->lgMeanAbundance() )
1004  continue;
1005 
1006  // loop over all ions of the isotope
1007  for( unsigned long ion=0; ion<it->second->ipMl.size(); ++ion )
1008  {
1009  if( it->second->ipMl[ion] != -1 &&
1010  (species[ it->second->ipMl[ion] ].location == NULL ) && (*atom)->ipMl[ion] != -1 )
1011  {
1012  species[ it->second->ipMl[ion] ].den = species[ (*atom)->ipMl[ion] ].den * it->second->frac;
1013  }
1014  }
1015  }
1016  }
1017 
1018  return;
1019 }
1020 
1022 {
1023  DEBUG_ENTRY( "t_mole_local::set_ion_locations()" );
1024 
1025  for( ChemNuclideList::iterator atom = nuclide_list.begin();
1026  atom != nuclide_list.end(); ++atom)
1027  {
1028  if ( !(*atom)->lgHasLinkedIon() )
1029  continue;
1030  long nelem = (*atom)->el->Z-1;
1031  for( long ion=0; ion<nelem+2; ++ion )
1032  {
1033  long mole_index = (*atom)->ipMl[ion];
1034  // element not enabled if index is -1
1035  if( mole_index != -1 )
1036  {
1037  ASSERT( mole_index < mole_global.num_total );
1038  species[mole_index].location = &(dense.xIonDense[nelem][ion]);
1039  dense.ionMole[nelem][ion] = &species[mole_index];
1040  }
1041  }
1042  }
1043 
1044  if( deut.lgElmtOn )
1045  {
1046  findspecieslocal("D")->location = &(deut.xIonDense[0]);
1047  findspecieslocal("D+")->location = &(deut.xIonDense[1]);
1048  }
1049 }
1050 
1052 {
1053  DEBUG_ENTRY( "total_molecule_deut()" );
1054 
1055  double total = 0.;
1056 
1057  if( !deut.lgElmtOn )
1058  return;
1059 
1060  for (long int i=0;i<mole_global.num_calc;++i)
1061  {
1062  if (mole.species[i].location == NULL && mole_global.list[i]->isIsotopicTotalSpecies() )
1063  {
1064  for( molecule::nNucsMap::iterator atom=mole_global.list[i]->nNuclide.begin();
1065  atom != mole_global.list[i]->nNuclide.end(); ++atom)
1066  {
1067  long int nelem = atom->first->el->Z-1;
1068  if( nelem==0 && atom->first->A==2 )
1069  {
1070  total += mole.species[i].den*atom->second;
1071  }
1072  }
1073  }
1074  }
1075 
1076  total_f = (realnum)total;
1077 
1078  return;
1079 }
1081 {
1082  DEBUG_ENTRY( "total_molecule_elems()" );
1083 
1084  /* now set total density of each element locked in molecular species */
1085  for ( long int nelem=ipHYDROGEN;nelem<LIMELM; ++nelem )
1086  {
1087  total[nelem] = 0.;
1088  }
1089  for (long int i=0;i<mole_global.num_calc;++i)
1090  {
1091  if (mole.species[i].location == NULL && mole_global.list[i]->isIsotopicTotalSpecies() )
1092  {
1093  for( molecule::nNucsMap::iterator atom=mole_global.list[i]->nNuclide.begin();
1094  atom != mole_global.list[i]->nNuclide.end(); ++atom)
1095  {
1096  ASSERT( atom->second > 0 );
1097  long int nelem = atom->first->el->Z-1;
1098  if( atom->first->lgMeanAbundance() )
1099  total[nelem] += (realnum) mole.species[i].den*atom->second;
1100  }
1101  }
1102  }
1103 }
1105 {
1106  long int i;
1107  realnum total;
1108 
1109  DEBUG_ENTRY( "total_molecules()" );
1110 
1111  total = 0.;
1112  for (i=0;i<mole_global.num_calc;++i)
1113  {
1114  if (mole.species[i].location == NULL && mole_global.list[i]->isIsotopicTotalSpecies())
1115  total += (realnum) mole.species[i].den;
1116  }
1117  return total;
1118 }
1120 {
1121  long int i;
1122  realnum total;
1123 
1124  DEBUG_ENTRY( "total_molecules_gasphase()" );
1125 
1126  total = 0.;
1127  for (i=0;i<mole_global.num_calc;++i)
1128  {
1129  if (mole_global.list[i]->lgGas_Phase && mole.species[i].location == NULL && mole_global.list[i]->isIsotopicTotalSpecies())
1130  total += (realnum) mole.species[i].den;
1131  }
1132  return total;
1133 }
1134 /*lint +e778 const express eval to 0 */
1135 /*lint +e725 expect positive indentation */
1136 
1138 {
1139  long int i, j;
1140  /* Neutrals and positive ions will be treated as single species inside
1141  molecular equilibrium solver, to facilitate coupling with ionization
1142  solvers */
1143  DEBUG_ENTRY ("mole_make_groups()" );
1144  if (mole_global.num_calc == 0)
1145  {
1146  groupspecies.resize(0);
1148  return;
1149  }
1151  for (i=0,j=0;i<mole_global.num_calc;i++)
1152  {
1153  if( mole_global.list[i]->isIsotopicTotalSpecies() && ( !mole_global.list[i]->isMonatomic() || mole_global.list[i]->charge <= 0 || ! mole_global.list[i]->lgGas_Phase ) )
1154  {
1155  /* Compound molecules and negative ions are represented individually */
1156  mole_global.list[i]->groupnum = j;
1157  groupspecies[j++] = &(*mole_global.list[i]);
1158  }
1159  else
1160  {
1161  /* All positive ions are collapsed into single macrostate (i.e. H+ -> H) */
1162  /* Need to increase constant if higher ions are included */
1163  ASSERT (mole_global.list[i]->charge < LIMELM+1);
1164  ASSERT (mole_global.list[i]->groupnum == -1);
1165  }
1166  }
1169 
1170  for (i=0;i<mole_global.num_calc;i++)
1171  {
1172  if (mole_global.list[i]->groupnum == -1)
1173  {
1174  if( mole_global.list[i]->isMonatomic() && mole_global.list[i]->isIsotopicTotalSpecies() )
1175  {
1176  for( nNucs_i it = mole_global.list[i]->nNuclide.begin(); it != mole_global.list[i]->nNuclide.end(); ++it )
1177  {
1178  ASSERT( it->second > 0 );
1179  mole_global.list[i]->groupnum = mole_global.list[it->first->ipMl[0]]->groupnum;
1180  break;
1181  }
1182  }
1183  else
1184  {
1185  ASSERT( !mole_global.list[i]->isIsotopicTotalSpecies() );
1186  mole_global.list[i]->groupnum = mole_global.list[ mole_global.list[i]->parentIndex ]->groupnum;
1187  }
1188  }
1189 
1190  ASSERT( mole_global.list[i]->groupnum != -1);
1191  }
1192 
1193  return;
1194 }
1195 
map< string, size_t >::iterator chem_nuclide_i
t_mole_global mole_global
Definition: mole.cpp:7
double grain_area
Definition: mole.h:457
nNucsMap nNuclide
Definition: mole.h:157
static void sort(MoleculeList::iterator start, MoleculeList::iterator end)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
const double ENERGY_H2_STAR
Definition: h2_priv.h:449
chem_element * el
Definition: mole.h:47
molecule::nNucsMap::iterator nNucs_i
Definition: mole.h:242
realnum getMass(int Anum)
Definition: abund.h:86
bool lgElmtOn
Definition: deuterium.h:21
molecule * null_mole
int n_nuclei(void) const
Definition: mole.h:161
vector< bool > lgTreatIsotopes
Definition: mole.h:359
int num_total
Definition: mole.h:362
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
mole_state
Definition: mole.h:19
realnum total_molecules(void)
int num_calc
Definition: mole.h:362
map< int, count_ptr< chem_nuclide > >::iterator isotopes_i
Definition: mole.h:39
void make_species(void)
void total_molecule_deut(realnum &total)
map< string, count_ptr< molecule > >::iterator molecule_i
Definition: mole_priv.h:40
bool lgLeiden_Keep_ipMH2s
Definition: hmi.h:219
bool lgTraceMole
Definition: trace.h:55
void updateXMolecules()
Definition: dense.cpp:24
const string label
Definition: mole.h:32
void mole_make_list(void)
map< string, count_ptr< mole_reaction > > reactab
STATIC molecule * newspecies(const char label[], enum spectype type, enum mole_state state, realnum form_enthalpy)
molezone * ionMole[LIMELM][LIMELM+1]
Definition: dense.h:137
bool speciesOff(const string &label)
FILE * ioQQQ
Definition: cddefines.cpp:7
vector< count_ptr< chem_element > > ChemElementList
Definition: mole.h:126
map< string, size_t > nuclidetab
molezone * findspecieslocal(const char buf[])
Definition: mole.h:378
bool lgGas_Phase
Definition: mole.h:160
void InitDeuteriumIonization()
Definition: deuterium.cpp:29
ChemNuclideList nuclide_list
Definition: mole.h:142
bool parse_species_label(const char label[], ChemNuclideList &atomsLeftToRight, vector< int > &numAtoms, string &embellishments)
int A
Definition: mole.h:48
int num_compacted
Definition: mole.h:362
realnum mass_amu
Definition: mole.h:51
molezone * findspecieslocal_validate(const char buf[])
t_dense dense
Definition: global.cpp:15
realnum form_enthalpy
Definition: mole.h:189
double grain_density
Definition: mole.h:457
t_elementnames elementnames
Definition: elementnames.cpp:5
Definition: mole.h:19
vector< int > ipMl
Definition: mole.h:50
ChemElementList element_list
bool isEnabled
Definition: mole.h:153
int getAiso(int iso)
Definition: abund.h:79
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
element_type * data() const
Definition: cddefines.h:1205
t_trace trace
Definition: trace.cpp:5
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:854
t_abund abund
Definition: abund.cpp:5
t_isotope IsoAbn[LIMELM]
Definition: abund.h:213
STATIC realnum MeanMassOfElement(const count_ptr< chem_element > &el)
molecule::nNucsMap::reverse_iterator nNucs_ri
Definition: mole.h:243
bool lgLeidenHack
Definition: mole.h:334
const ios_base::openmode mode_r
Definition: cpu.h:267
map< int, count_ptr< chem_nuclide > > isotopes
Definition: mole.h:33
STATIC bool isactive(const molecule &mol)
#define STATIC
Definition: cddefines.h:118
map< string, count_ptr< chem_element > > elemtab
string parentLabel
Definition: mole.h:147
STATIC void read_species_file(string filename, bool lgCreateIsotopologues=true)
int groupnum
Definition: mole.h:194
chem_nuclide * null_nuclide
t_mole_local mole
Definition: mole.cpp:8
molecule * findspecies(const char buf[])
void create_isotopologues_one_position(unsigned position, ChemNuclideList &atoms, vector< int > &numAtoms, string atom_old, string atom_new, string embellishments, string &newLabel)
t_atoms atoms
Definition: atoms.cpp:5
int getNiso()
Definition: abund.h:66
float realnum
Definition: cddefines.h:124
valarray< class molezone > species
Definition: mole.h:468
#define EXIT_FAILURE
Definition: cddefines.h:168
enum mole_state state
Definition: mole.h:193
count_ptr< chem_nuclide > findnuclide(const char buf[])
STATIC void SetIsotopeFractions(const vector< bool > &lgResolveNelem)
#define cdEXIT(FAIL)
Definition: cddefines.h:482
int index
Definition: mole.h:194
realnum total_molecules_gasphase(void)
map< string, count_ptr< mole_reaction > > functab
const int Z
Definition: mole.h:31
STATIC void newelement(const char label[], int ipion)
realnum mole_return_cached_species(const GroupMap &MoleMap)
bool lgExcit
Definition: mole.h:159
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
int n_react
Definition: mole.h:195
bool exists(const molecule *m)
Definition: mole.h:258
void SetDeuteriumFractionation(const realnum &frac)
Definition: deuterium.cpp:57
void SetGasPhaseDeuterium(const realnum &Hdensity)
Definition: deuterium.cpp:65
t_state state
Definition: state.cpp:18
double xIonDense[2]
Definition: deuterium.h:23
bool lgElmtOn[LIMELM]
Definition: dense.h:160
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
int compare(const molecule &mol2) const
Definition: mole.h:211
realnum gas_phase[LIMELM]
Definition: dense.h:76
vector< count_ptr< chem_nuclide > > ChemNuclideList
Definition: mole.h:124
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
#define ASSERT(exp)
Definition: cddefines.h:613
T * get_ptr() const
Definition: count_ptr.h:49
void set_isotope_abundances(void)
STATIC void start(long, realnum[], realnum[], long, long[], realnum *, int *)
realnum getAbn(int Anum)
Definition: abund.h:121
Definition: abund.h:39
const int LIMELM
Definition: cddefines.h:308
bool isMonatomic(void) const
Definition: mole.h:182
int charge
Definition: mole.h:158
double den
Definition: mole.h:421
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
STATIC void newisotope(const count_ptr< chem_element > &el, int massNumberA, realnum mass_amu, double frac)
bool lgGrain_mole_deplete
Definition: mole.h:356
bool lgDifferByExcitation(const molecule &mol1, const molecule &mol2)
vector< GrainBin * > bin
Definition: grainvar.h:583
void mole_update_species_cache(void)
STATIC bool ispassive(const molecule &mol)
t_deuterium deut
Definition: deuterium.cpp:7
string label
Definition: mole.h:156
double elec
Definition: mole.h:460
realnum mole_mass
Definition: mole.h:190
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
MoleculeList list
Definition: mole.h:365
vector< molecule * > groupspecies
double frac
Definition: mole.h:52
spectype
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
const int ipCARBON
Definition: cddefines.h:354
string label(void) const
Definition: mole.h:66
bool lgNoHeavyMole
Definition: mole.h:328
#define fixit(a)
Definition: cddefines.h:417
bool lgElemsConserved(void)
Definition: dense.cpp:119
GrainVar gv
Definition: grainvar.cpp:5
t_hmi hmi
Definition: hmi.cpp:5
double frac(double d)
void set_ion_locations()
bool lgNoMole
Definition: mole.h:325
bool lgDustOn() const
Definition: grainvar.h:475
const int ipHYDROGEN
Definition: cddefines.h:349
void mole_make_groups(void)
map< string, count_ptr< molecule > > spectab
chem_element * null_element
void total_molecule_elems(realnum total[LIMELM])
molecule * findspecies_validate(const char buf[])
molezone * null_molezone
double * location
Definition: mole.h:411
void create_isotopologues(ChemNuclideList &atoms, vector< int > &numAtoms, string atom_old, string atom_new, string embellishments, vector< string > &newLabels)
double grain_saturation
Definition: mole.h:457
void updateXMolecules()
Definition: deuterium.cpp:16