Hi,
I will be more specific, starting the problem since begining:
1)I have these images (map of force due a potential, Fx and Fy) that I need to replicate the last two , using mathematica ( my adviser wanted to use Fortran..., but I am fighting to consolidate Mathematica coding in my department )
The potential that I need to Minimize at each step is:
-V *(2 Cos[(2 Pi x)/a]*
Cos[(2 Pi y)/(a Sqrt[3])] + Cos[(4 Pi y)/(a Sqrt[3])]) + (k/2)*(X- x)^2 + (k/2)*(Y - y)^2
So, when I get the minima, or (x,y) positions, they will be used to the next step. I use them like that:
first step:
i(=X) = 0 and j (=Y)=0, -> potential minimized -> x,y calculated.
second step:
i= 0 j=1, potential minimized using the last (x,y).....-> "new" x,y calculated..
and so on...
And I am going collecting at each step the pairs (x,y) , that will be used to calculate the Map Forces Fx and Fy.
At the final I have an output.txt with all x,y, collected. and I use x for calculate Fx and y for calculate Fy.
and finally I make a table to plot Fx x X(or i) and Fy x Y(or j)
obs: "F = k* x" (oscilator)
well, and I did a lot of tests (since august....) and the best images I could get (using contourplot) were (for Fx, and Fy, respectively)
The parameters:
(*a= 2.46 \[Angstrom] = 2.46 10^-10m*)
(*k= 25 N/m as 1 joule = N.m, so , N/m= eV/(16 \[Angstrom]^2)*)
And the code is that:
$HistoryLength = 0;
ClearSystemCache[];
ClearAll["Global`*"];
SetDirectory["/home/marcelo/Documentos/Tomlinson/"];
Block[{x, y, i, j, FmolaX, FmolaY, c = {0, 0}, a = 2.46, k = 1.5625,
V = 0.5},
Emin[go_, ho_] :=
Module[{g = go, h = ho}, {x, y} /.
Take[Flatten[
N[FindMinimum[-V *(2 Cos[(2 Pi x)/a]*
Cos[(2 Pi y)/(a Sqrt[3])] +
Cos[(4 Pi y)/(a Sqrt[3])]) + (k/2)*(x - g)^2 + (k/
2)*(y - h)^2, {{x, First[c]}, {y, Last[c]}},
Method -> "PrincipalAxis"], 7]], -2]];
Do[
If[j == 0, c = {0, i}]; (*<- to avoid contour problems in the beginning of the step, I can´t use the last x,y results, for the next i*)
c = Emin[j, i] ;
FmolaX = -16*k*(First[c] - (j))*(10^-10);
FmolaY = -16*k*(Last[c] - (i))*(10^-10);
{FmolaX, FmolaY} >>> "DataFxFyTomlinsonV20.txt";
c >>> "DataMinimTomlinsonV20.txt",
(*Print[
"the minima of the function related to step Y",i," X",j,":",
c],*) {i, 0.0, 10.0, 0.1}, {j, 0.0, 10.0, 0.1}]
];
If you can see, I did steps of 0.1, and each step I used the last pair of {x,y},of the minimized potential.
As my images were not so good, mainly Fy, I tried aplying the gradient :
F=-\[Del](x,y) Subscript[E, tot]
So , I got two eqautions from the potential, and then I tested NSolve for Fx=0 and Fy = 0.
Details, When I tested the model 1d worked fine! See below:
equauni =
Quiet[NSolve[(2 \[Pi] )/ a Subscript[V, z] Sin[(2 \[Pi] x)/a] +
k (x - #) == 0, x, Reals] & /@
Table[i, {i, 0, 10, 0.02}] /. {a -> 246/100, k -> 25/16,
Subscript[V, z] -> 5/10}] >> datateste1
equauni2 = Flatten[Take[{x} /. equauni // Simplify, All, 1]];
dataequani = Transpose[{Table[i, {i, 0, 10, 0.02}], equauni2}];
Map[Length, {equauni2}];
Map[Length, {Table[i, {i, 0, 10, 0.02}]}];
ListPlot[Tooltip[Transpose[{Table[i, {i, 0, 10, 0.02}], equauni2}]],
PlotStyle -> Green,
AxesLabel -> {Subscript[x, 0][\[Angstrom]], x[\[Angstrom]]}]
And calculating the Force:
Fmola[x_,
d_] := -k (x -
d)*(10^-10);(*converted to Newton (N)*)
ListPlot[Tooltip[
Transpose[{ 10^-10 equauni2,
Fmola[equauni2, Table[i, {i, 0, 10, 0.02}]] /. {a -> 246/100,
k -> 25/1, Subscript[V, z] -> 5/10}}]],
AxesLabel -> {x["m"], F["N"]}, AxesOrigin -> {0, 0},
PlotRange -> All, ImageSize -> Large, PlotStyle -> Green, Mesh -> All]
I hope, this explanation can be more useful.
