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 package ffx.algorithms.commands.test;
39
40 import ffx.algorithms.cli.AlgorithmsCommand;
41 import ffx.numerics.Potential;
42 import ffx.potential.bonded.Atom;
43 import ffx.potential.bonded.Molecule;
44 import ffx.utilities.FFXBinding;
45 import picocli.CommandLine.Command;
46 import picocli.CommandLine.Option;
47 import picocli.CommandLine.Parameters;
48
49 import java.util.ArrayList;
50 import java.util.Arrays;
51 import java.util.List;
52
53 import static java.lang.String.format;
54
55
56
57
58
59
60
61
62
63 @Command(description = " Find guest atoms to restrain near host molecule.", name = "test.FindRestraints")
64 public class FindRestraints extends AlgorithmsCommand {
65
66
67
68
69 @Option(names = {"--hostName"}, paramLabel = "BCD", defaultValue = "BCD",
70 description = "Host molecule name in the file.")
71 private String hostName;
72
73
74
75
76 @Option(names = {"--guestName"}, paramLabel = "LIG", defaultValue = "LIG",
77 description = "Ligand molecule name in the file.")
78 private String guestName;
79
80
81
82
83 @Option(names = {"--distanceCutoff"}, paramLabel = "5", defaultValue = "5",
84 description = "Cutoff to use when selecting guest atoms near host COM")
85 private double distanceCutoff;
86
87
88
89
90 @Parameters(arity = "1..*", paramLabel = "files",
91 description = "XYZ or PDB input files.")
92 private List<String> filenames;
93
94
95
96
97
98 public Potential potential;
99
100
101
102
103
104 public Potential getPotentialObject() {
105 return potential;
106 }
107
108
109
110
111 public FindRestraints() {
112 super();
113 }
114
115
116
117
118
119 public FindRestraints(FFXBinding binding) {
120 super(binding);
121 }
122
123
124
125
126
127 public FindRestraints(String[] args) {
128 super(args);
129 }
130
131
132
133
134 @Override
135 public FindRestraints run() {
136 if (!init()) {
137 return this;
138 }
139
140 activeAssembly = getActiveAssembly(filenames.get(0));
141 if (activeAssembly == null) {
142 logger.info(helpString());
143 return this;
144 }
145
146 Molecule[] molArr = activeAssembly.getMoleculeArray();
147
148 List<Atom> restrainHostList = new ArrayList<>();
149 List<Atom> restrainList = new ArrayList<>();
150 double[] COM = new double[3];
151 double[] subCOM = new double[3];
152 int[] restrainHostIndices = new int[]{11, 16, 17, 20, 23, 26, 31, 32, 39, 40, 51, 63, 64, 70, 71, 82, 94, 95, 101, 102, 113, 125, 126, 132, 133, 144, 156, 157, 163, 164, 175, 187, 188, 191, 198};
153 for (Molecule molecule : molArr) {
154 logger.info(format(" Molecule name: " + molecule.getName()));
155 if (molecule.getName().contains(hostName)) {
156 Atom[] host_atoms = molecule.getAtomList().toArray(new Atom[0]);
157 COM = getCOM(host_atoms);
158 logger.info(format(" Center of mass of host molecule: " + Arrays.toString(COM)));
159 for (Atom atom : host_atoms) {
160 if (contains(restrainHostIndices, atom.getIndex())) {
161 restrainHostList.add(atom);
162 logger.info(format(" Atom: " + atom));
163 }
164 }
165 Atom[] subAtoms = restrainHostList.toArray(new Atom[0]);
166 subCOM = getCOM(subAtoms);
167 logger.info(format(" Center of mass of subsection host atoms: " + Arrays.toString(subCOM)));
168 double comdist = Math.sqrt(Math.pow(subCOM[0] - COM[0], 2) +
169 Math.pow(subCOM[1] - COM[1], 2) +
170 Math.pow(subCOM[2] - COM[2], 2));
171 logger.info(format(" Distance between COMs: " + comdist));
172 } else if (molecule.getName().contains(guestName)) {
173 Atom[] guest_atoms = molecule.getAtomList().toArray(new Atom[0]);
174 for (Atom atom : guest_atoms) {
175 double dist = Math.sqrt(Math.pow(atom.getXYZ().get()[0] - subCOM[0], 2) +
176 Math.pow(atom.getXYZ().get()[1] - subCOM[1], 2) +
177 Math.pow(atom.getXYZ().get()[2] - subCOM[2], 2));
178
179
180
181
182 if (dist < distanceCutoff && atom.isHeavy()) {
183 restrainList.add(atom);
184 }
185 }
186 }
187 }
188 logger.info(format(" Number of atoms to restrain: " + restrainList.size()));
189 int[] restrainIndices = restrainList.stream()
190 .mapToInt(Atom::getIndex)
191 .toArray();
192 logger.info(format(" Restrain list indices: " + Arrays.toString(restrainIndices)));
193
194 return this;
195 }
196
197
198
199
200
201
202
203 private static double[] getCOM(Atom[] atoms) {
204
205 double[] COM = new double[3];
206 double totalMass = 0.0;
207 for (Atom s : atoms) {
208 double[] pos = s.getXYZ().get();
209 COM[0] += pos[0] * s.getMass();
210 COM[1] += pos[1] * s.getMass();
211 COM[2] += pos[2] * s.getMass();
212 totalMass += s.getMass();
213 }
214 totalMass = 1 / totalMass;
215 COM[0] *= totalMass;
216 COM[1] *= totalMass;
217 COM[2] *= totalMass;
218
219 return COM;
220 }
221
222
223
224
225
226
227
228
229 private static boolean contains(int[] array, int value) {
230 return Arrays.stream(array).anyMatch(element -> element == value);
231 }
232 }