mardi 1 octobre 2019

Can someone Interpret python code to the c++

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