So for Part 1 I simply had to make two vectors X & Y and make them so that when

Question

So for Part 1 I simply had to make two vectors X & Y and make them so that when plotted ( scatter(x,y);) they were each 3.8 units apart in a plot and now im not sure how to do part 2. Thanks

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) 4& i 1 [606 120 or 12 LU 4& xi 13 ij where r V(xi xi)2 lyi yj)2, or 3.361 Angstrom 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 3vanderWaals (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.

Explanation / Answer

function [energy,Fx,Fy] = vanderWalls(x,y)
   clc; clear all; close all
   VR = linspace(x,y,100);
   figure(1);
   ylim([0 2])
   xlim([x y])
   xlabel('VR')
   ylabel('PR')
   Tr = 0.9;
   Prfunc = @(VR) 8*Tr./(3*VR - 1) - 3./(VR.^2);
   PR = Prfunc(VR);
   plot(VR,PR)
   if Tr < 1
       energy = 1.0;
       VDW_res = [1 -1/3*(1+8*Tr/energy) 3/energy -1/energy];
       res = sort(roots(VDW_res));
       Fx = (res(2)-res(1))*energy - integral(Prfunc,res(1),res(2));
       Fy = integral(Prfunc,res(2),res(3)) - (res(3)-res(2))*energy;
       Z = abs(Fx-Fy);
       while Z > 0.0001
           VDW_res = [1 -1/3*(1+8*Tr/energy) 3/energy -1/energy];
           res = sort(roots(VDW_res));
           Prfunc = @(VR) 8*Tr./(3*VR - 1) - 3./(VR.^2);
           Fx = (res(2)-res(1))*energy - integral(Prfunc,res(1),res(2));
           Fy = integral(Prfunc,res(2),res(3)) - (res(3)-res(2))*energy;
           Z = abs(Fx - Fy);
           energy = energy - 0.00001;
           figure(1); hold off;
           plot(VR,PR)
           figure(1); hold on;
           plot([x y],[energy energy],'k--')
           hold off;
       end
   end
end
[energy,Fx,Fy] = vanderWalls(0.25,3);
disp(energy)
disp(Fx)
disp(Fy)

Related Questions

Navigate

• Browse (All)
• Browse S
• Subjects

---

---

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