From f73b0abaed63851cb7b1b0bd659d3469f2576691 Mon Sep 17 00:00:00 2001 From: Ondrej Marsalek Date: Thu, 8 Jan 2015 11:22:27 -0800 Subject: [PATCH 1/2] Add PIMD NM mode - all NM frequencies are omega_n --- ipi/engine/normalmodes.py | 34 ++++++++++++++++++++++++---------- ipi/inputs/normalmodes.py | 6 +++--- 2 files changed, 27 insertions(+), 13 deletions(-) diff --git a/ipi/engine/normalmodes.py b/ipi/engine/normalmodes.py index 38566736..fc77c1c9 100644 --- a/ipi/engine/normalmodes.py +++ b/ipi/engine/normalmodes.py @@ -273,26 +273,28 @@ def get_nmm(self): # also checks that the frequencies and the mode given in init are # consistent with the beads and ensemble - dmf = np.zeros(self.nbeads,float) - dmf[:] = 1.0 + sk = np.ones(self.nbeads) + if self.mode == "rpmd": if len(self.nm_freqs) > 0: - warning("nm.frequencies will be ignored for RPMD mode.", verbosity.low) + warning("Normal mode frequencies will be ignored for RPMD mode.", verbosity.low) elif self.mode == "manual": if len(self.nm_freqs) != self.nbeads-1: raise ValueError("Manual path mode requires (nbeads-1) frequencies, one for each internal mode of the path.") for b in range(1, self.nbeads): - sk = self.omegak[b]/self.nm_freqs[b-1] - dmf[b] = sk**2 + sk[b] = self.omegak[b]/self.nm_freqs[b-1] + elif self.mode == "pimd": + if len(self.nm_freqs) > 0: + warning("Normal mode frequencies will be ignored for PIMD mode.", verbosity.low) + for b in range(1, self.nbeads): + sk[b] = self.omegak[b] / self.omegan elif self.mode == "pa-cmd": if len(self.nm_freqs) > 1: warning("Only the first element in nm.frequencies will be considered for PA-CMD mode.", verbosity.low) if len(self.nm_freqs) == 0: raise ValueError("PA-CMD mode requires the target frequency of all the internal modes.") for b in range(1, self.nbeads): - sk = self.omegak[b]/self.nm_freqs[0] - info(" ".join(["NM FACTOR", str(b), str(sk), str(self.omegak[b]), str(self.nm_freqs[0])]), verbosity.medium) - dmf[b] = sk**2 + sk[b] = self.omegak[b]/self.nm_freqs[0] elif self.mode == "wmax-cmd": if len(self.nm_freqs) > 2: warning("Only the first two element in nm.frequencies will be considered for WMAX-CMD mode.", verbosity.low) @@ -301,8 +303,20 @@ def get_nmm(self): wmax = self.nm_freqs[0] wt = self.nm_freqs[1] for b in range(1, self.nbeads): - sk = 1.0/np.sqrt((wt)**2*(1+(wmax/self.omegak[1])**2)/(wmax**2+(self.omegak[b])**2)) - dmf[b] = sk**2 + sk[b] = 1.0/np.sqrt((wt)**2*(1+(wmax/self.omegak[1])**2)/(wmax**2+(self.omegak[b])**2)) + + dmf = sk**2 + + # TODO + # - nicely print info at medium verbosity level regardless of mode + # - do this elsewhere so that it is not repeated - aparently it gets printed twice here + # - also print the actual frequency used in dynamics - is that just `get_dynwk` above? + # - replace original code: + # info(" ".join(["NM FACTOR", str(b), str(sk), str(self.omegak[b]), str(self.nm_freqs[0])]), verbosity.medium) + print 'DBG | NM | mode =', self.mode + print 'DBG | NM | omegan = ', self.omegan + for b in range(self.nbeads): + print 'DBG | NM | %3d %12.6f %12.6f %12.6f' % (b, sk[b], dmf[b], self.omegak[b]) return dmf diff --git a/ipi/inputs/normalmodes.py b/ipi/inputs/normalmodes.py index 0e460510..13c95f64 100644 --- a/ipi/inputs/normalmodes.py +++ b/ipi/inputs/normalmodes.py @@ -9,7 +9,7 @@ This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License @@ -43,8 +43,8 @@ class InputNormalModes(InputArray): attribs = copy(InputArray.attribs) attribs["mode"] = (InputAttribute, {"dtype" : str, "default" : "rpmd", - "help" : "Specifies the technique to be used to calculate the dynamical masses. 'rpmd' simply assigns the bead masses the physical mass. 'manual' sets all the normal mode frequencies except the centroid normal mode manually. 'pa-cmd' takes an argument giving the frequency to set all the non-centroid normal modes to. 'wmax-cmd' is similar to 'pa-cmd', except instead of taking one argument it takes two ([wmax,wtarget]). The lowest-lying normal mode will be set to wtarget for a free particle, and all the normal modes will coincide at frequency wmax. ", - "options" : ['pa-cmd', 'wmax-cmd', 'manual', 'rpmd']}) + "help" : "Specifies the technique to be used to calculate the dynamical masses. 'rpmd' simply assigns the bead masses the physical mass. 'pimd' sets all the normal mode frequencies to omega_n. 'manual' sets all the normal mode frequencies except the centroid normal mode manually. 'pa-cmd' takes an argument giving the frequency to set all the non-centroid normal modes to. 'wmax-cmd' is similar to 'pa-cmd', except instead of taking one argument it takes two ([wmax,wtarget]). The lowest-lying normal mode will be set to wtarget for a free particle, and all the normal modes will coincide at frequency wmax. ", + "options" : ['pa-cmd', 'wmax-cmd', 'manual', 'rpmd', 'pimd']}) attribs["transform"] = (InputValue,{"dtype" : str, "default" : "fft", "help" : "Specifies whether to calculate the normal mode transform using a fast Fourier transform or a matrix multiplication. For small numbers of beads the matrix multiplication may be faster.", From 4b2c0a0ce9b0751b787d3a38a73c85cd9eae8ac6 Mon Sep 17 00:00:00 2001 From: Ondrej Marsalek Date: Fri, 9 Jan 2015 12:56:47 -0800 Subject: [PATCH 2/2] Correct target frequency; text output still hackish As the dynamics runs at nT, the potential should also be scaled down by n to achieve the same effect as in PIMD with temperature T and target frequency omega_P. --- ipi/engine/normalmodes.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/ipi/engine/normalmodes.py b/ipi/engine/normalmodes.py index fc77c1c9..b2fe76e2 100644 --- a/ipi/engine/normalmodes.py +++ b/ipi/engine/normalmodes.py @@ -287,7 +287,7 @@ def get_nmm(self): if len(self.nm_freqs) > 0: warning("Normal mode frequencies will be ignored for PIMD mode.", verbosity.low) for b in range(1, self.nbeads): - sk[b] = self.omegak[b] / self.omegan + sk[b] = self.omegak[b] / (self.omegan / self.nbeads) elif self.mode == "pa-cmd": if len(self.nm_freqs) > 1: warning("Only the first element in nm.frequencies will be considered for PA-CMD mode.", verbosity.low) @@ -307,18 +307,24 @@ def get_nmm(self): dmf = sk**2 + return dmf + + def get_nm_report(self): + """ """ + # TODO + # - return, rather than print + # - add a header # - nicely print info at medium verbosity level regardless of mode # - do this elsewhere so that it is not repeated - aparently it gets printed twice here - # - also print the actual frequency used in dynamics - is that just `get_dynwk` above? # - replace original code: # info(" ".join(["NM FACTOR", str(b), str(sk), str(self.omegak[b]), str(self.nm_freqs[0])]), verbosity.medium) print 'DBG | NM | mode =', self.mode - print 'DBG | NM | omegan = ', self.omegan + print 'DBG | NM | omega_n = ', self.omegan + for i, f in enumerate(self.nm_freqs): + print 'DBG | NM | %3d %12.6f' % (i ,f) for b in range(self.nbeads): - print 'DBG | NM | %3d %12.6f %12.6f %12.6f' % (b, sk[b], dmf[b], self.omegak[b]) - - return dmf + print 'DBG | NM | %3d %12.6f %12.6f %12.6f' % (b, self.omegak[b], self.nm_factor[b], self.dynomegak[b]) def get_dynm3(self): """Returns an array with the dynamical masses of individual atoms in the normal modes representation.""" @@ -327,6 +333,9 @@ def get_dynm3(self): for b in range(self.nbeads): dm3[b] = self.beads.m3[b]*self.nm_factor[b] + # TODO: call this from somewhere else! + self.get_nm_report() + return dm3 def free_qstep(self):