-
Notifications
You must be signed in to change notification settings - Fork 0
/
read.go
153 lines (128 loc) · 3.37 KB
/
read.go
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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
package gr
import (
"bufio"
"fmt"
"io"
"strconv"
"strings"
"github.com/kpotier/molsolvent/pkg/util"
)
// readCfgFirst reads the first configuration. It reads the number of atoms, the
// columns and performs the usual calculations like in readCfg.
func (g *GR) readCfgFirst(r *bufio.Reader) (box [3]float64, xyz XYZ, err error) {
g.atoms, box, err = util.Header(r, nil, readSlice)
if err != nil {
err = fmt.Errorf("Header: %w", err)
return
}
b, _ := r.ReadSlice('\n')
fields := strings.Fields(string(b))
if len(fields) <= 2 {
err = fmt.Errorf("not enough columns (at least 3; got %d)", len(fields))
return
}
fields = fields[2:]
var found int
g.colsLen = len(fields)
for k, v := range fields {
switch v {
case "x":
g.cols[0] = k
case "y":
g.cols[1] = k
case "z":
g.cols[2] = k
case "type":
g.cols[3] = k
default:
continue
}
found++
}
if found < len(g.cols) {
return box, nil, fmt.Errorf("cannot find the columns x, y, z, and type")
}
g.order, xyz, err = g.fetchXYZFirst(r)
if err != nil {
return box, nil, fmt.Errorf("fetchXYZ: %w", err)
}
return
}
// readCfg reads a configuration of the LAMMPS trajectory. This method will call
// fetchXYZ to fetch the coordinates of the two atoms.
func (g *GR) readCfg(r *bufio.Reader) (box [3]float64, xyz XYZ, err error) {
box, err = util.HeaderWOutAtoms(r, nil, readSlice)
if err != nil {
err = fmt.Errorf("HeaderWOutAtoms: %w", err)
return
}
r.ReadSlice('\n')
xyz, err = g.fetchXYZ(r)
if err != nil {
err = fmt.Errorf("fetchXYZ: %w", err)
return
}
return
}
// fetchXYZ fetches the coordinates of the two atoms by calling readXYZ two
// times (one for the first atom, and the other for the second atom). This
// method is like fetchXYZ but it returns the order of the atoms.
func (g *GR) fetchXYZFirst(r *bufio.Reader) (order []string, xyz XYZ, err error) {
xyz = make(XYZ, len(g.atomsTyp))
nbat := g.atoms / len(g.atomsTyp)
for _, v := range g.atomsTyp {
xyz[v] = make([][3]float64, 0, nbat)
}
for i := 0; i < g.atoms; i++ {
var typ string
typ, err = g.readXYZ(r, xyz)
if err != nil {
return
}
if _, ok := g.Atoms[typ]; ok {
order = append(order, typ)
}
}
return
}
// fetchXYZ fetches the coordinates of the two atoms by calling readXYZ two
// times (one for the first atom, and the other for the second atom).
func (g *GR) fetchXYZ(r *bufio.Reader) (xyz XYZ, err error) {
xyz = make(map[string][][3]float64, len(g.atomsTyp))
nbat := g.atoms / len(g.atomsTyp)
for _, v := range g.atomsTyp {
xyz[v] = make([][3]float64, 0, nbat)
}
for i := 0; i < g.atoms; i++ {
_, err = g.readXYZ(r, xyz)
if err != nil {
return
}
}
return
}
// readXYZ reads the coordinates for each atom. If the atom type exists in XYZ,
// it is added to the map. It returns the type of the atom.
func (g *GR) readXYZ(r *bufio.Reader, xyz XYZ) (typ string, err error) {
b, _ := r.ReadSlice('\n')
fields := strings.Fields(string(b))
if len(fields) != g.colsLen {
err = fmt.Errorf("number of columns don't match: %d (expected %d)", len(fields), g.colsLen)
return
}
typ = fields[g.cols[3]]
xyzTyp, ok := xyz[typ]
if !ok {
return
}
var xyzTmp [3]float64
for k := 0; k < 3; k++ {
xyzTmp[k], _ = strconv.ParseFloat(fields[g.cols[k]], 64)
}
xyz[typ] = append(xyzTyp, xyzTmp)
return
}
func readSlice(r *bufio.Reader, w io.Writer) []byte {
b, _ := r.ReadSlice('\n')
return b
}