Simpatico  v1.10
ChainMaker.cpp
1 #include "ChainMaker.h"
2 #include <util/format/Int.h>
3 #include <util/format/Dbl.h>
4 
5 #include <string>
6 
7 namespace Tools
8 {
9 
10  ChainMaker::ChainMaker(){};
11 
12  ChainMaker::~ChainMaker(){};
13 
14  void ChainMaker::readParam(std::istream& in)
15  {
16  bondPotential_.setNBondType(1);
17 
18  read<Boundary>(in, "boundary", boundary_);
19  readParamComposite(in, bondPotential_);
20  readParamComposite(in, random_);
21  read<int>(in, "nMolecule", nMolecule_);
22  read<int>(in, "nAtomPerMolecule", nAtomPerMolecule_);
23 
24  in >> outputStyle_;
25  }
26 
27 
28  void ChainMaker::writeChains(std::ostream& out)
29  {
30  if (outputStyle_ == "McMd") {
31  writeChainsMcMd(out);
32  } else
33  if (outputStyle_ == "DdMd") {
34  writeChainsDdMd(out);
35  } else
36  if (outputStyle_ == "DdMdMolecule") {
37  writeChainsDdMdMole(out);
38  } else {
39  std::cout << "Unrecognized style "
40  << outputStyle_ << std::endl;
41  }
42  }
43 
44  void ChainMaker::writeChainsMcMd(std::ostream& out)
45  {
46  Vector r;
47  Vector v;
48  double beta = 1.0;
49  int bondType = 0;
50  int iMol, iAtom;
51 
52  out << "BOUNDARY" << std::endl;
53  out << std::endl;
54  out << boundary_ << std::endl;
55  out << std::endl;
56  out << "MOLECULES" << std::endl;
57  out << std::endl;
58  out << "species " << 0 << std::endl;
59  out << "nMolecule " << nMolecule_ << std::endl;
60  out << std::endl;
61  for (iMol = 0; iMol < nMolecule_; ++iMol) {
62  out << "molecule " << iMol << std::endl;
63  // Choosen first atom at random
64  boundary_.randomPosition(random_, r);
65  out << r << std::endl;
66  for (iAtom = 1; iAtom < nAtomPerMolecule_; ++iAtom) {
67  random_.unitVector(v);
68  v *= bondPotential_.randomBondLength(&random_, beta, bondType);
69  r += v;
70  boundary_.shift(r);
71  out << r << std::endl;
72  }
73  out << std::endl;
74  }
75 
76  }
77 
78  void ChainMaker::writeChainsDdMd(std::ostream& out)
79  {
80  Vector r;
81  Vector velocity;
82  Vector v;
83  double beta = 1.0;
84  int atomType = 0;
85  int bondType = 0;
86  int iMol, iAtom, i, j;
87 
88  out << "BOUNDARY" << std::endl;
89  out << std::endl;
90  out << boundary_ << std::endl;
91  out << std::endl;
92  out << "ATOMS" << std::endl;
93  out << "nAtom " << nMolecule_*nAtomPerMolecule_ << std::endl;
94 
95  i = 0;
96  velocity.zero();
97  for (iMol = 0; iMol < nMolecule_; ++iMol) {
98 
99  boundary_.randomPosition(random_, r);
100 
101  for (iAtom = 0; iAtom < nAtomPerMolecule_; ++iAtom) {
102 
103  out << Int(i,6) << Int(atomType, 10);
104  for (j = 0; j < Dimension; ++j) {
105  out << Dbl(r[j], 15, 6);
106  }
107  out << " ";
108  for (j = 0; j < Dimension; ++j) {
109  out << Dbl(0.0, 12, 4);
110  }
111  out << std::endl;
112 
113  if (iAtom < nAtomPerMolecule_ - 1) {
114  random_.unitVector(v);
115  v *= bondPotential_.randomBondLength(&random_, beta, bondType);
116  r += v;
117  boundary_.shift(r);
118  }
119 
120  ++i;
121  }
122  }
123 
124  // Write bonds
125  out << std::endl;
126  out << "BONDS" << std::endl;
127  out << "nBond " << nMolecule_*(nAtomPerMolecule_ -1 ) << std::endl;
128  i = 0;
129  j = 0;
130  for (iMol = 0; iMol < nMolecule_; ++iMol) {
131  for (iAtom = 0; iAtom < nAtomPerMolecule_ - 1; ++iAtom) {
132  out << Int(j,5) << Int(bondType, 5) << " ";
133  out << Int(i, 10) << Int(i + 1, 10) << std::endl;
134  ++i;
135  ++j;
136  }
137  ++i;
138  }
139 
140  // Write angles
141  int angleType = 0;
142  out << std::endl;
143  out << "ANGLES" << std::endl;
144  out << "nAngle " << nMolecule_*(nAtomPerMolecule_ -2) << std::endl;
145  i = 0;
146  j = 0;
147  for (iMol = 0; iMol < nMolecule_; ++iMol) {
148  for (iAtom = 0; iAtom < nAtomPerMolecule_ - 2; ++iAtom) {
149  out << Int(j,5) << Int(angleType, 5) << " ";
150  out << Int(i, 10) << Int(i + 1, 10) << Int(i + 2, 10)
151  << std::endl;
152  ++i;
153  ++j;
154  }
155  i += 2;
156  }
157 
158  // Write dihedrals
159  int dihedralType = 0;
160  out << std::endl;
161  out << "DIHEDRALS" << std::endl;
162  out << "nDihedral " << nMolecule_*(nAtomPerMolecule_ -3) << std::endl;
163  i = 0;
164  j = 0;
165  for (iMol = 0; iMol < nMolecule_; ++iMol) {
166  for (iAtom = 0; iAtom < nAtomPerMolecule_ - 3; ++iAtom) {
167  out << Int(j,5) << Int(dihedralType, 5) << " ";
168  out << Int(i, 10) << Int(i + 1, 10) << Int(i + 2, 10)
169  << Int(i + 3, 10) << std::endl;
170  ++i;
171  ++j;
172  }
173  i += 3;
174  }
175 
176  }
177 
178  void ChainMaker::writeChainsDdMdMole(std::ostream& out)
179  {
180  Vector r;
181  Vector velocity;
182  Vector v;
183  double beta = 1.0;
184  int atomType = 0;
185  int sId = 0;
186  int bondType = 0;
187  int iMol, iAtom, i, j;
188 
189  out << "BOUNDARY" << std::endl;
190  out << std::endl;
191  out << boundary_ << std::endl;
192  out << std::endl;
193  out << "ATOMS" << std::endl;
194  out << "nAtom " << nMolecule_*nAtomPerMolecule_ << std::endl;
195 
196  i = 0;
197  velocity.zero();
198  for (iMol = 0; iMol < nMolecule_; ++iMol) {
199 
200  boundary_.randomPosition(random_, r);
201 
202  for (iAtom = 0; iAtom < nAtomPerMolecule_; ++iAtom) {
203 
204  out << Int(i,10) << Int(atomType, 6)
205  << Int(sId,6) << Int(iMol, 10) << Int(iAtom,6);
206  out << "\n";
207  for (j = 0; j < Dimension; ++j) {
208  out << Dbl(r[j], 15, 6);
209  }
210  out << "\n";
211  for (j = 0; j < Dimension; ++j) {
212  out << Dbl(0.0, 15, 6);
213  }
214  out << "\n";
215 
216  if (iAtom < nAtomPerMolecule_ - 1) {
217  random_.unitVector(v);
218  v *= bondPotential_.randomBondLength(&random_, beta, bondType);
219  r += v;
220  boundary_.shift(r);
221  }
222 
223  ++i;
224  }
225  out << std::endl;
226  }
227 
228  // Write bonds
229  out << std::endl;
230  out << "BONDS" << std::endl;
231  out << "nBond " << nMolecule_*(nAtomPerMolecule_ -1 ) << std::endl;
232  i = 0;
233  j = 0;
234  for (iMol = 0; iMol < nMolecule_; ++iMol) {
235  for (iAtom = 0; iAtom < nAtomPerMolecule_ - 1; ++iAtom) {
236  out << Int(j,5) << Int(bondType, 5) << " ";
237  out << Int(i, 10) << Int(i + 1, 10) << std::endl;
238  ++i;
239  ++j;
240  }
241  ++i;
242  }
243 
244  // Write angles
245  int angleType = 0;
246  out << std::endl;
247  out << "ANGLES" << std::endl;
248  out << "nAngle " << nMolecule_*(nAtomPerMolecule_ -2) << std::endl;
249  i = 0;
250  j = 0;
251  for (iMol = 0; iMol < nMolecule_; ++iMol) {
252  for (iAtom = 0; iAtom < nAtomPerMolecule_ - 2; ++iAtom) {
253  out << Int(j,5) << Int(angleType, 5) << " ";
254  out << Int(i, 10) << Int(i + 1, 10) << Int(i + 2, 10)
255  << std::endl;
256  ++i;
257  ++j;
258  }
259  i += 2;
260  }
261 
262  // Write dihedrals
263  int dihedralType = 0;
264  out << std::endl;
265  out << "DIHEDRALS" << std::endl;
266  out << "nDihedral " << nMolecule_*(nAtomPerMolecule_ -3) << std::endl;
267  i = 0;
268  j = 0;
269  for (iMol = 0; iMol < nMolecule_; ++iMol) {
270  for (iAtom = 0; iAtom < nAtomPerMolecule_ - 3; ++iAtom) {
271  out << Int(j,5) << Int(dihedralType, 5) << " ";
272  out << Int(i, 10) << Int(i + 1, 10) << Int(i + 2, 10)
273  << Int(i + 3, 10) << std::endl;
274  ++i;
275  ++j;
276  }
277  i += 3;
278  }
279 
280  }
281 
282 }
283 
284 int main()
285 {
286  Tools::ChainMaker obj;
287  obj.readParam(std::cin);
288  obj.writeChains(std::cout);
289 }
Util::Vector::zero
Vector & zero()
Set all elements of a 3D vector to zero.
Definition: Vector.h:514
Util::Dimension
const int Dimension
Dimensionality of space.
Definition: Dimension.h:19
Util::Vector
A Vector is a Cartesian vector.
Definition: Vector.h:75
Simp::OrthorhombicBoundary::randomPosition
void randomPosition(Random &random, Vector &r) const
Generate random position within the primary unit cell.
Definition: OrthorhombicBoundary.cpp:105
Tools::ChainMaker::readParam
void readParam(std::istream &in)
Read the parameter file block.
Definition: ChainMaker.cpp:14
Util::Dbl
Wrapper for a double precision number, for formatted ostream output.
Definition: Dbl.h:39
Tools::ChainMaker
Generates a simulation box full of randomly generated linear polymers.
Definition: ChainMaker.h:21
Simp::HarmonicBond::setNBondType
void setNBondType(int nBondType)
Set the number of bond types.
Definition: HarmonicBond.cpp:57
Util::Int
Wrapper for an int, for formatted ostream output.
Definition: Int.h:36
Simp::OrthorhombicBoundary::shift
void shift(Vector &r) const
Shift Cartesian Vector r to its primary image.
Definition: OrthorhombicBoundary.h:432
Tools
Single-processor classes for pre- and post-processing MD trajectories.
Definition: tools/analyzers/Analyzer.cpp:13
Simp::HarmonicBond::randomBondLength
double randomBondLength(Random *random, double beta, int type) const
Return bond length chosen from equilibrium distribution.
Definition: HarmonicBond.cpp:120
Util::ParamComposite::readParamComposite
void readParamComposite(std::istream &in, ParamComposite &child, bool next=true)
Add and read a required child ParamComposite.
Definition: ParamComposite.cpp:260
Util::Random::unitVector
void unitVector(Vector &v)
Generate unit vector with uniform probability over the unit sphere.
Definition: Random.cpp:122
Generated on Wed Mar 7 2018 12:59:38 for Simpatico by   doxygen 1.8.11