md.py

import matplotlib.pyplot as plt
import random as rd

import initial as init, verlet as verl
from properties import *
from input import *
from metadynamics import *
from processing import *

# Main Loop.
# ------------------------------------------------------------------------

def simulate():
    # -----------------------Initialize--------------------------------
    ## system
    N = Ncube ** 3
    R = init.InitPositionCubic(Ncube, L)
    V = init.InitVelocity(N, T0, M)
    E = np.zeros(steps)

    ## metadynamics
    n_gauss = 0
    S = [] # position of the center of the Gaussian

    print("steps\ttemperature\tpressure\tkinetic\tpotential\tenergy\t")

    for t in range(0, steps):
        # -----------------------Measuring Physical Properties----------------------------

        ## calculate kinetic energy contribution
        k = my_kinetic_energy(V, M)

        ## calculate temperature
        T = my_temperature(k, N)

        ## calculate distance table
        drij = get_distance_table(N, R, L)

        ## calculate potential energy contribution
        p = my_potential_energy(drij, rc)

        ## calculate total energy
        E[t] = k + p

        ## calculate forces
        F = np.array([my_force_on(i, R, L, rc) for i in range(N)])
        A = F / M

        ## calculate pressure
        P = my_pressure(N**3,N,T,R,F)

        # -----------------------Anderson Thermostat----------------------
        if anderson == True:
            sigma = (Ta / M) ** 0.5
            mean = 0
            for i in V:
                if np.random.random() < eta * h:
                    i[0] = rd.gauss(mean, sigma)
                    i[1] = rd.gauss(mean, sigma)
                    i[2] = rd.gauss(mean, sigma)

        # -----------------------Propagation----------------------------
        ## calculate new positions
        nR = verl.VerletNextR(R, V, A, h)
        nR = my_pos_in_box(nR, L)

        ## calculate forces with new positions nR
        nF = np.array([my_force_on(i, nR, L, rc) for i in range(N)])

        ## bias forces with metadynamics
        metaF = meta(t, n_gauss, S, nF, nR)

        ## calculate new velocities
        nA = nF / M
        nV = verl.VerletNextV(V, A, nA, h)

        # update positions:
        R, V = nR, nV

        # ------------------------Output-------------------------------------
        print('%d\t%.3f\t%.2e\t%.5f\t%.5f\t%.5f\t' % (t, T, P, k, p, E[t]))

    return S


if __name__ == '__main__':
    # Run the simulation
    S = simulate()

    # Post processing
    postprocessing(S)