function [energy, fx, fy] = vanderWaals(x, y) eps4 = 0.2824 * 4.0; rmin = 3.361;

Question

function [energy, fx, fy] = vanderWaals(x, y)
eps4 = 0.2824 * 4.0;
rmin = 3.361;
rmin6 = rmin^6;
rmin12 = rmin^12;
nAtoms = length(x);
% Zero out energy.
energy = 0;
% Zero out forces.
fx = zeros(1,nAtoms);
fy = zeros(1,nAtoms);
% Loop over atoms to accumulate energy and forces.
for i=1:nAtoms
xi = x(i);
yi = y(i);
% Loop over atoms whose index is greater than i.
for j=i+1:nAtoms
% Compute separation distance (Angstroms).
dx = xi - x(j);
dy = yi - y(j);
r = sqrt(dx*dx + dy*dy);
% Compute interaction energy.
ir = 1.0/r;
ir6 = ir^6;
ir12 = ir6*ir6;
e = eps4*(rmin12*ir12 - rmin6*ir6);
energy = energy + e;
% Compute equal and opposite forces, includeing chain rule term.
ir7 = ir6*ir;
ir13 = ir12*ir;
de = eps4*(6.0*rmin6*ir7 - 12.0*rmin12*ir13);
dxr = de*dx/r;
dyr = de*dy/r;
fx(i) = fx(i) - dxr;
fy(i) = fy(i) - dyr;
fx(j) = fx(j) + dxr;
fy(j) = fy(j) + dyr;
end
end
end

Part II: Initial Ener Forces and Velocities. Use the supplied Matlab function vanderWaals to compute the initial potential energy and forces on each atom of the Argon system. This approximates atomic interactions by summing the van der Waals energy between all pairs of atoms as: 12 ECX) j i +1 i 1 606 120 12 or 4E xi 13 ax where rij V (xi xi)2 (yi yi)2, o 3.361 Angstroms and E 0.2824 Kcal/mole. The two inputs are the coordinate vectors 'x' and 'y' and the returned values are the energy (kcal/mole) and forces (kcal/mole/A): [energy, FX, Fyl vanderWaals (x,y); Print out the returned energy. Finally, compute the initial velocities in the x and y directions using "randn"with mean 0.0 and std. dev. 0.05. HINT: The energy of the Argon cluster for the initial Cartesian coordinates should be exactly -13.8298 Kcal/mole. Part III: Simulation of the Argon Cluster for 5000 Time Steps 1. From the positions, velocities and atomic forces of the ith step determine the coordinates for step (i+1). Assume m-1 so that Fx max simplifies to Fx ax and choose a time step At 0.01 (psec): xi vriAt ax i At 2/2 Xi+1 yi vyiAt ay,i At 2/2 yi+1 2. Using the new coordinates, call the vanderWaals function to update the potential energy and forces. Finally, update the velocities: vri (ari x, i+1 x, i+1 vyi (ay i At/2 vy,i +1 y,i+1 Note: The new velocities depend on current and previous accelerations. 3. Compute the kinetic energy (kcal/mole) as: kin 1 4. Repeat Steps 1-3 for 5000 time steps. Save the potential energy and kinetic energy at each time step (their sum is constant due to conservation of energy)

Explanation / Answer

function [energy, fx, fy] = vanderWaals(x, y)
eps4 = 0.2824 * 4.0;
rmin = 3.361;
rmin6 = rmin^6;
rmin12 = rmin^12;
nAtoms = length(x);
% Zero out energy.
energy = 0;
% Zero out forces.
fx = zeros(1,nAtoms);
fy = zeros(1,nAtoms);
% Loop over atoms to accumulate energy and forces.
for i=1:nAtoms
xi = x(i);
yi = y(i);
% Loop over atoms whose index is greater than i.
for j=i+1:nAtoms
% Compute separation distance (Angstroms).
dx = xi - x(j);
dy = yi - y(j);
r = sqrt(dx*dx + dy*dy);
% Compute interaction energy.
ir = 1.0/r;
ir6 = ir^6;
ir12 = ir6*ir6;
e = eps4*(rmin12*ir12 - rmin6*ir6);
energy = energy + e;
% Compute equal and opposite forces, includeing chain rule term.
ir7 = ir6*ir;
ir13 = ir12*ir;
de = eps4*(6.0*rmin6*ir7 - 12.0*rmin12*ir13);
dxr = de*dx/r;
dyr = de*dy/r;
fx(i) = fx(i) - dxr;
fy(i) = fy(i) - dyr;
fx(j) = fx(j) + dxr;
fy(j) = fy(j) + dyr;
end
end
end

Related Questions

Navigate

• Browse (All)
• Browse F
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.