-
Notifications
You must be signed in to change notification settings - Fork 0
/
ext.cpp
117 lines (91 loc) · 3.09 KB
/
ext.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#include "cellFile.h"
#include "input.h"
#include "ext.h"
void Extend::Routine()
{
TITLE("Extend","Routine");
// (1) Initialization
// cel_in : input geometry file
// cel_out: output geometry file
CellFile cel_in;
CellFile cel_out;
cel_in.file_name = INPUT.geo_in;
// (2) Read in the geometry
if( !CellFile::ReadGeometry( cel_in ) )
{
cout << " Can't find the geometry file : " << INPUT.geo_in << endl;
exit(0);
}
// (3) Extend the cell.
Extend::ExtendCell( cel_in, cel_out );
bool cartesian = false;
cartesian = INPUT.cartesian;
// make sure this value is set.
cel_in.ntype = INPUT.ntype;
cel_out.ntype = cel_in.ntype;
// (4) Output the geometry.
cout << " should write geometry here" << endl;
CellFile::WriteGeometry( cel_out, cartesian );
return;
}
void Extend::ExtendCell( const Cell &cel_in, Cell &cel_out )
{
TITLE("Extend","ExtendCell");
// (1) initialization
cel_out.system_name = cel_in.system_name;
cel_out.coordinate = cel_in.coordinate;
cout << " cel_in.system_name=" << cel_in.system_name << endl;
const int ntype = INPUT.ntype;
cel_out.atom = new Atoms[ntype];
// (2) expand the lattice vectors
cel_out.a1.x = cel_in.a1.x * INPUT.ext_1;
cel_out.a1.y = cel_in.a1.y * INPUT.ext_1;
cel_out.a1.z = cel_in.a1.z * INPUT.ext_1;
cel_out.a2.x = cel_in.a2.x * INPUT.ext_2;
cel_out.a2.y = cel_in.a2.y * INPUT.ext_2;
cel_out.a2.z = cel_in.a2.z * INPUT.ext_2;
cel_out.a3.x = cel_in.a3.x * INPUT.ext_3;
cel_out.a3.y = cel_in.a3.y * INPUT.ext_3;
cel_out.a3.z = cel_in.a3.z * INPUT.ext_3;
// (3) calculate the output atom positions.
Vector3<double> add1,add2,add3;
bool frac = true;
int ia2=0;
cel_out.nat = 0;
cout << " Extension = " << INPUT.ext_1 << " * " << INPUT.ext_2 << " * " << INPUT.ext_3 << endl;
for(int it=0; it<ntype; ++it)
{
cel_out.atom[it].id = cel_in.atom[it].id;
cel_out.atom[it].pot_file = cel_in.atom[it].pot_file;
// this sentence has problem
cel_out.atom[it].na = cel_in.atom[it].na * INPUT.ext_1 * INPUT.ext_2 * INPUT.ext_3;
cout << " number of atoms in new cell = " << cel_out.atom[it].na << endl;
cel_out.nat += cel_out.atom[it].na;
cel_out.atom[it].posd = new Vector3<double>[cel_out.atom[it].na];
cel_out.atom[it].pos = new Vector3<double>[cel_out.atom[it].na];
ia2 = 0; // mohan fix bug 2013-06-24
for(int ia=0; ia<cel_in.atom[it].na; ++ia)
{
for(int i=0; i<INPUT.ext_1; i++)
{
for(int j=0; j<INPUT.ext_2; j++)
{
for(int k=0; k<INPUT.ext_3; k++)
{
// calculate the new atom positions using fractional coordinates.
cel_out.atom[it].posd[ia2].x = cel_in.atom[it].posd[ia].x + (double)i;
cel_out.atom[it].posd[ia2].y = cel_in.atom[it].posd[ia].y + (double)j;
cel_out.atom[it].posd[ia2].z = cel_in.atom[it].posd[ia].z + (double)k;
cel_out.atom[it].posd[ia2].x /= (double)INPUT.ext_1;
cel_out.atom[it].posd[ia2].y /= (double)INPUT.ext_2;
cel_out.atom[it].posd[ia2].z /= (double)INPUT.ext_3;
cel_out.direct2cartesian(it,ia2);
++ia2;
}
}
}
}
}
cout << " After extension, the total atom number is " << cel_out.nat << endl;
return;
}