I want to interpret python code to the C++, to improve the execution speed. The code is about polymer physics and Monte carlo methods. I also wrote a bunch of C++ expressions related to python expressions but having the difficulty to put them in order. Here is the python code:
import numpy as np
def gen_chain(N):
coordinates = np.loadtxt('saw.txt', skiprows=0)
return coordinates
def lj(rij2):
sig_by_r6 = np.power(sigma / rij2, 3)
sig_by_r12 = np.power(sig_by_r6, 2)
lje = 4.0 * epsilon * (sig_by_r12 - sig_by_r6)
return lje
def fene(rij2):
return (-0.5 * K * R**2 * np.log(1 - ((np.sqrt(rij2) - r0)**2 / R**2)))
def total_energy(coord):
# Non-bonded
e_nb = 0
for i in range(N):
for j in range(i - 1):
ri = coord[i]
rj = coord[j]
rij = ri - rj
rij2 = np.dot(rij, rij)
if (np.sqrt(rij2) < rcutoff):
e_nb += lj(rij2)
# Bonded
e_bond = 0
for i in range(1, N):
ri = coord[i]
rj = coord[i - 1]
rij = ri - rj
rij2 = np.dot(rij, rij)
e_bond += fene(rij2)
return e_nb + e_bond
def move(coord):
trial = np.ndarray.copy(coord)
for i in range(N):
delta = (2.0 * np.random.rand(3) - 1) * max_delta
trial[i] += delta
return trial
def accept(delta_e):
beta = 1.0 / T
if delta_e <= 0.0:
return True
random_number = np.random.rand(1)
p_acc = np.exp(-beta * delta_e)
if random_number < p_acc:
return True
return False
if __name__ == "__main__":
# FENE parameters
K = 40
R = 0.3
r0 = 0.7
# LJ parameters
sigma = 0.624
epsilon = 1.0
# MC parameters
N = 100 # number of particles
rcutoff = 2.5 * sigma
max_delta = 0.01
n_steps = 100000
T = 0.5
coord = gen_chain(N)
energy_current = total_energy(coord)
traj = open('traj.xyz', 'w')
for step in range(n_steps):
if step % 1000 == 0:
traj.write(str(N) + '\n\n')
for i in range(N):
traj.write("C %10.5f %10.5f %10.5f\n" % (coord[i][0], coord[i][1], coord[i][2]))
print(step, energy_current)
coord_trial = move(coord)
energy_trial = total_energy(coord_trial)
delta_e = energy_trial - energy_current
if accept(delta_e):
coord = coord_trial
energy_current = energy_trial
traj.close()
Below you will find the C++ code:
double kb = 1;
double sigma = 1;
double epsilon = 1;
double dt; // discrete time step
double t;
int N;
void accept(const double t){ // Aceept or reject
double r, Enew;
Enew=energy();
r=ran->doub();
if (r < exp(-(Enew-E)/t) ) {// accept
xyz=xyzn;
E=Enew;
}
return;
}
double LJ_calc (double x, double xs, double ep){
double x_rel = xs / x;
double LJP = ep *( pow( x_rel , 12 ) - 2 * pow( x_rel , 6 ) );
return LJP;
}
// A Finite Extendible Non-linear Elastic potential
double U_fene(double r2, double k, double a){
if( r2<a*a ){
return -0.5*k*a*a*log(1.0-r2/a/a);
}
}
I don't want to put the Monte carlo step because it is confusing. I hope someone can answer. I know it will be very time consuming.
Aucun commentaire:
Enregistrer un commentaire