Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Post equilibrium analysis #43

Open
wants to merge 15 commits into
base: main
Choose a base branch
from
Open

Post equilibrium analysis #43

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
778d917
commiting so I can pull changes
Mar 29, 2023
b6860ac
Added rdf code to project.py. Still needs testing.
Apr 12, 2023
7f6d5d9
rdf code working
Apr 17, 2023
0623171
added code for average PE
Apr 18, 2023
2939e12
commiting so I can push changes
Apr 18, 2023
0a3be67
Merge branch 'post-eq' of github.com:cmelab/Advanced-Sampling into po…
Apr 18, 2023
eff0f0d
fixing some errors I left
Apr 18, 2023
b74534e
added code that calculates the mean PE and adds it to the job documen…
Apr 20, 2023
8af4ee9
Merge branch 'post-eq' of github.com:cmelab/Advanced-Sampling into po…
Apr 20, 2023
25cc2fd
removing useless files from this branch of the repo
Apr 20, 2023
ef79380
Deleting Workspace
gwenwhite Apr 20, 2023
45cbee4
removing workspace directory so changes can be pulled
Apr 20, 2023
c141637
added average nearest neighbor calculation. Also tidied up some things.
Apr 20, 2023
2094721
chaning the rmax value to the default value
Apr 20, 2023
efe2073
added standard deviation of PE calculation
Apr 20, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 5 additions & 4 deletions MCMC/MCMC-flow/init.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,13 +34,13 @@ def get_parameters():
parameters["m"] = [6]

# logging parameters
parameters["energy_write_freq"] = [1000]
parameters["energy_write_freq"] = [20000]
parameters["trajectory_write_freq"] = [10000]

# run parameters
parameters["n_steps"] = [[1e7, 1e8]]
parameters["kT"] = [[10, 1.5]]
parameters["max_trans"] = [[3.0, 0.5]]
parameters["n_steps"] = [[1e7, 2e8]]
parameters["kT"] = [[10, 1.0]]
parameters["max_trans"] = [[3.0, 0.3]]
parameters["seed"] = [20]

return list(parameters.keys()), list(product(*parameters.values()))
Expand All @@ -59,6 +59,7 @@ def main():
parent_job.init()
parent_job.doc.setdefault("done", False)
parent_job.doc.setdefault("timestep", [])
parent_job.doc.setdefault("analyzed", False)
parent_job.doc.setdefault("accepted_moves", [])
parent_job.doc.setdefault("rejected_moves", [])
parent_job.doc.setdefault("acceptance_ratio", [])
Expand Down
30 changes: 30 additions & 0 deletions MCMC/MCMC-flow/project.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,9 @@ def sampled(job):
def initialized(job):
return job.isfile("trajectory.gsd")

@MyProject.label
def analyzed(job):
return job.doc.get("analyzed")

@directives(executable="python -u")
@MyProject.operation
Expand Down Expand Up @@ -131,6 +134,33 @@ def sample(job):
print("Simulation finished completed")
print("-----------------------------")

@MyProject.operation
@MyProject.pre(sampled)
@MyProject.post(analyzed)
def analysis(job):
from utils import avg_nn, rdf
import numpy as np
os.makedirs(os.path.join(job.ws, "analysis/rdf"))
gsdfile = job.fn('trajectory_1.gsd')
rdf = rdf(gsdfile, start=-30)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please apply this change to the rdf function in utils.py

rdf.compute((box,points), reset=False)

The default value for reset is True, which means if you are computing rdf from multiple snapshots, freud doesn't accumulate the computated values from different snapshots. When reset is True, the result rdf is only from last snapshot.

x = rdf.bin_centers
y = rdf.rdf
peakx = max(x)
peaky = max(y)
Comment on lines +148 to +149
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think first and second peaks are more information rich in rdf. The max of bin_centers is always going to be the last value in the list (largest distance) and max of rdf doesn't give us a correct estimation of the first peak.
Previously I used scipy.signal.find_peaks from scipy package to find the peak.
@chrisjonesBSU Feel free to correct me if I'm wrong.

Copy link
Member

@chrisjonesBSU chrisjonesBSU May 9, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree. We should use some kind of peak finding algorithm or tool. RDFs often have more than 1 peak, so just finding the max y value may not provide enough information.

logfile = job.fn('log.txt')
energy = np.genfromtxt(logfile, delimiter=",")
mean = np.mean(energy[:,0])
sd = np.std(energy[:,0])
save_path = os.path.join(job.ws, "analysis/rdf/rdf.txt")
np.savetxt(save_path, np.transpose([x,y]), delimiter=',', header ="bin_centers, rdf")
save_peak = os.path.join(job.ws, "analysis/rdf/peak.txt")
np.savetxt(save_peak, np.transpose([peakx, peaky]), delimiter=',', header="max_x, max_y")
job.doc['avg_PE'] = mean
job.doc['sd_PE'] = sd
nn = avg_nn(gsdfile, frame=-1)
job.doc['average_nn'] = nn
job.doc["analyzed"] = True


if __name__ == "__main__":
MyProject().main()
10 changes: 10 additions & 0 deletions MCMC/log.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
344.14088824424385,1.0
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

let's add .txt to the .gitignore and remove this file.

344.14088824424385,1.0
341.59388377626556,1.0
338.6685686445551,1.0
337.94650985867054,1.0
337.9376561616747,1.0
336.81233332171166,1.0
333.6200195542917,1.0
332.86849037229774,1.0
331.1718868176446,1.0