1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
9 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
10 //
11 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
12 //////////////////////////////////////////////////////////////////////////////////////
13
14
15 #include "ParticleIOUtility.h"
16 #include "Utilities/ProgressReportEngine.h"
17 namespace qmcplusplus
18 {
19 #if OHMMS_DIM == 3
expandSuperCell(ParticleSet & ref_,Tensor<int,3> & tmat)20 void expandSuperCell(ParticleSet& ref_, Tensor<int, 3>& tmat)
21 {
22 typedef ParticleSet::SingleParticlePos_t SingleParticlePos_t;
23 typedef ParticleSet::Tensor_t Tensor_t;
24 Tensor<int, 3> I(1, 0, 0, 0, 1, 0, 0, 0, 1);
25 bool identity = true;
26 int ij = 0;
27 while (identity && ij < 9)
28 {
29 identity = (I[ij] == tmat[ij]);
30 ++ij;
31 }
32 if (identity)
33 return;
34 ReportEngine PRE("expandSuperCell", " ");
35 app_log() << " TileMatrix != Identity. Expanding a simulation cell for " << ref_.getName() << std::endl;
36 {
37 char buff[500];
38 snprintf(buff, 500, " tilematrix= %4d %4d %4d %4d %4d %4d %4d %4d %4d\n", tmat[0], tmat[1], tmat[2], tmat[3],
39 tmat[4], tmat[5], tmat[6], tmat[7], tmat[8]);
40 app_log() << buff << std::endl;
41 }
42 //convert2unit
43 ref_.convert2Unit(ref_.R);
44 ParticleSet::ParticleLayout_t PrimCell(ref_.Lattice);
45 ref_.Lattice.set(dot(tmat, PrimCell.R));
46 int natoms = ref_.getTotalNum();
47 int numCopies = std::abs(det(tmat));
48 ParticleSet::ParticlePos_t primPos(ref_.R);
49 ParticleSet::ParticleIndex_t primTypes(ref_.GroupID);
50 ref_.resize(natoms * numCopies);
51 int maxCopies = 10;
52 int index = 0;
53 //set the unit to the Cartesian
54 ref_.R.InUnit = PosUnit::Cartesian;
55 app_log() << " Reduced coord Cartesion coord species.\n";
56 for (int ns = 0; ns < ref_.getSpeciesSet().getTotalNum(); ++ns)
57 {
58 for (int i0 = -maxCopies; i0 <= maxCopies; i0++)
59 for (int i1 = -maxCopies; i1 <= maxCopies; i1++)
60 for (int i2 = -maxCopies; i2 <= maxCopies; i2++)
61 for (int iat = 0; iat < primPos.size(); iat++)
62 {
63 if (primTypes[iat] != ns)
64 continue;
65 //SingleParticlePos_t r = primPos[iat];
66 SingleParticlePos_t uPrim = primPos[iat];
67 for (int i = 0; i < 3; i++)
68 uPrim[i] -= std::floor(uPrim[i]);
69 SingleParticlePos_t r = PrimCell.toCart(uPrim) + (double)i0 * PrimCell.a(0) + (double)i1 * PrimCell.a(1) +
70 (double)i2 * PrimCell.a(2);
71 SingleParticlePos_t uSuper = ref_.Lattice.toUnit(r);
72 if ((uSuper[0] >= -1.0e-6) && (uSuper[0] < 0.9999) && (uSuper[1] >= -1.0e-6) && (uSuper[1] < 0.9999) &&
73 (uSuper[2] >= -1.0e-6) && (uSuper[2] < 0.9999))
74 {
75 char buff[500];
76 snprintf(buff, 500, " %10.4f %10.4f %10.4f %12.6f %12.6f %12.6f %d\n", uSuper[0], uSuper[1],
77 uSuper[2], r[0], r[1], r[2], ns);
78 app_log() << buff;
79 ref_.R[index] = r;
80 ref_.GroupID[index] = ns; //primTypes[iat];
81 ref_.ID[index] = index;
82 ref_.PCID[index] = iat;
83 index++;
84 }
85 }
86 }
87 app_log() << " Simulationcell after tiling" << std::endl;
88 ref_.Lattice.print(app_log());
89 app_log() << std::endl;
90 }
91
92 #elif OHMMS_DIM == 2
93 void expandSuperCell(ParticleSet& ref_, Tensor<int, 2>& tmat)
94 {
95 typedef ParticleSet::SingleParticlePos_t SingleParticlePos_t;
96 typedef ParticleSet::Tensor_t Tensor_t;
97 Tensor<int, 2> I(1, 0, 0, 1);
98 bool identity = true;
99 int ij = 0;
100 while (identity && ij < 4)
101 {
102 identity = (I[ij] == tmat[ij]);
103 ++ij;
104 }
105 if (identity)
106 return;
107 //convert2unit
108 ref_.convert2Unit(ref_.R);
109 ParticleSet::ParticleLayout_t PrimCell(ref_.Lattice);
110 ref_.Lattice.set(dot(tmat, PrimCell.R));
111 int natoms = ref_.getTotalNum();
112 int numCopies = std::abs(det(tmat));
113 ParticleSet::ParticlePos_t primPos(ref_.R);
114 ParticleSet::ParticleIndex_t primTypes(ref_.GroupID);
115 ref_.resize(natoms * numCopies);
116 int maxCopies = 10;
117 int index = 0;
118 //set the unit to the Cartesian
119 ref_.R.InUnit = PosUnit::CartesianUnit;
120 app_log() << " Reduced coord Cartesion coord species.\n";
121 for (int ns = 0; ns < ref_.getSpeciesSet().getTotalNum(); ++ns)
122 {
123 for (int i0 = -maxCopies; i0 <= maxCopies; i0++)
124 for (int i1 = -maxCopies; i1 <= maxCopies; i1++)
125 for (int iat = 0; iat < primPos.size(); iat++)
126 {
127 if (primTypes[iat] != ns)
128 continue;
129 //SingleParticlePos_t r = primPos[iat];
130 SingleParticlePos_t uPrim = primPos[iat];
131 for (int i = 0; i < 2; i++)
132 uPrim[i] -= std::floor(uPrim[i]);
133 SingleParticlePos_t r = PrimCell.toCart(uPrim) + (double)i0 * PrimCell.a(0) + (double)i1 * PrimCell.a(1);
134 SingleParticlePos_t uSuper = ref_.Lattice.toUnit(r);
135 if ((uSuper[0] >= -1.0e-6) && (uSuper[0] < 0.9999) && (uSuper[1] >= -1.0e-6) && (uSuper[1] < 0.9999))
136 {
137 char buff[500];
138 snprintf(buff, 500, " %10.4f %10.4f %12.6f %12.6f %d\n", uSuper[0], uSuper[1], r[0], r[1], ns);
139 app_log() << buff;
140 ref_.R[index] = r;
141 ref_.GroupID[index] = ns; //primTypes[iat];
142 index++;
143 }
144 }
145 }
146 }
147 #else
148 #error "Only 2D and 3D are implemented in ParticleIO/ParticleIOUtilitcy.cpp"
149 #endif
150 } // namespace qmcplusplus
151