Commit 9727fdc4 authored by julient31's avatar julient31
Browse files

Commit JT 110818 

- correct bug (match ewald/disp results for vir)
- started correct mag. part
parent d5fe8857

examples/SPIN/pppm_spin/data.2

0 → 100644
+13 −0
Original line number Diff line number Diff line 
RANDOM INITIALIZATION FOR STOCKMAYER FLUID

2     atoms

1     atom types



     -3.0   3.0  xlo xhi

     -3.0   3.0  ylo yhi     

     -3.0   3.0  zlo zhi

     #30.0 30.0 0.0    xy xz yz



Atoms



1	1	1.73	 0.0	 0.0	 0.0	 0.0	 0.0	 1.0 

2	1	1.73	 0.0	 2.5	 0.0	 0.0	 0.0	 1.0

examples/SPIN/pppm_spin/in.spin.2

0 → 100644
+75 −0
Original line number Diff line number Diff line 
# two magnetic atoms in a 3d box



clear 

units	 	metal

atom_style 	spin



dimension 	3

#boundary 	p p p

atom_modify 	map array 



read_data	../examples/SPIN/pppm_spin/data.2





mass		1 58.93

#set 		group all spin/random 31 1.72



#velocity 	all create 100 4928459 rot yes dist gaussian



pair_style 	spin/dipolar/cut 4.0

pair_coeff 	* * long 2.6 

#pair_style 	hybrid/overlay spin/exchange 4.0 spin/dipolar/long 8.0

#pair_style 	hybrid/overlay eam/alloy spin/exchange 4.0 spin/dipolar/long 8.0

#pair_style 	hybrid/overlay eam/alloy spin/exchange 4.0 spin/dipolar/long/qsymp 8.0

#pair_coeff 	* * eam/alloy ../examples/SPIN/pppm_spin/Co_PurjaPun_2012.eam.alloy Co

#pair_coeff 	* * spin/exchange exchange 4.0 0.1 1.135028015e-05 1.064568567

#pair_coeff 	* * spin/dipolar/long long 8.0



#neighbor 	0.1 bin

#neigh_modify 	every 10 check yes delay 20

neighbor	0.3 bin 

neigh_modify	delay 0

#neigh_modify	every 1 delay 10 check yes page 100000000 one 10000000



#kspace_style 	pppm/dipole/spin 1.0e-4

#kspace_style 	ewald/dipole/spin 1.0e-4

#kspace_modify 	compute yes

#kspace_modify 	compute yes gewald 0.1



fix 		1 all precession/spin zeeman 0.0 0.0 0.0 1.0

fix 		2 all langevin/spin 0.0 0.0 21

#fix 		3 all nve/spin lattice yes

fix 		3 all nve/spin lattice no



timestep	0.0001



thermo_style	custom step temp pe ke etotal press

thermo_modify   format float %20.16g

thermo		1



#compute		peratom all pe/atom

#compute		pe all reduce sum c_peratom

#thermo_style	custom step temp pe c_pe



#compute peratom2 all stress/atom

#compute peratom2 all stress/atom NULL

#compute p all reduce sum c_peratom2[1] c_peratom2[2] c_peratom2[3]  c_peratom2[4] c_peratom2[5] c_peratom2[6] 

#variable press equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol)

#variable pxx equal -c_p[1]/vol

#variable pyy equal -c_p[2]/vol

#variable pzz equal -c_p[3]/vol

#variable pxy equal -c_p[4]/vol

#variable pxz equal -c_p[5]/vol

#variable pyz equal -c_p[6]/vol

#thermo_style custom step temp etotal pe c_pe press v_press pxx v_pxx pyy v_pyy pzz v_pzz pxy v_pxy pxz v_pxz pyz v_pyz

#thermo_style custom step etotal pe press v_press v_pxx v_pyy v_pzz v_pxy v_pxz v_pyz

#thermo_style custom step temp etotal press v_press



compute 	outsp all property/atom spx spy spz sp fmx fmy fmz

dump		1 all custom 1 dump.equil id type x y z c_outsp[1] c_outsp[2] c_outsp[3] 

#c_outsp[5] c_outsp[6] c_outsp[7]

#dump_modify 1 format line "%d %d %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g" scale yes



#pair_modify compute no



run	 	1

examples/SPIN/pppm_spin/in.spin.cut_comp

0 → 100644
+72 −0
Original line number Diff line number Diff line 
# bcc iron in a 3d periodic box



clear 

units	 	metal

atom_style 	spin



dimension 	3

boundary 	p p p

atom_modify 	map array 



lattice 	bcc 2.8665

region 		box block 0.0 5.0 0.0 5.0 0.0 5.0

create_box 	1 box

create_atoms 	1 box



mass		1 58.93

#set 		group all spin 2.2 0.0 0.0 1.0

set 		group all spin/random 31 2.2



#pair_style 	spin/dipolar/cut 5.0

#pair_coeff 	* * long 5.0

pair_style 	spin/dipolar/long 4.0

pair_coeff 	* * long 4.0



#neighbor 	0.1 bin

#neigh_modify 	every 10 check yes delay 20

neighbor	0.3 bin 

neigh_modify	delay 0

#neigh_modify	every 1 delay 10 check yes page 100000000 one 10000000



#kspace_style 	pppm/dipole/spin 1.0e-4

kspace_style 	ewald/dipole/spin 1.0e-4

kspace_modify 	compute yes

#kspace_modify 	compute yes gewald 0.1



fix 		1 all precession/spin zeeman 0.0 0.0 0.0 1.0

fix 		2 all langevin/spin 0.0 0.1 21

#fix 		3 all nve/spin lattice yes

fix 		3 all nve/spin lattice no



timestep	0.0001



thermo_style	custom step temp pe ke etotal press

thermo_modify   format float %20.16g

thermo		50



#compute		peratom all pe/atom

#compute		pe all reduce sum c_peratom

#thermo_style	custom step temp pe c_pe



#compute peratom2 all stress/atom

#compute peratom2 all stress/atom NULL

#compute p all reduce sum c_peratom2[1] c_peratom2[2] c_peratom2[3]  c_peratom2[4] c_peratom2[5] c_peratom2[6] 

#variable press equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol)

#variable pxx equal -c_p[1]/vol

#variable pyy equal -c_p[2]/vol

#variable pzz equal -c_p[3]/vol

#variable pxy equal -c_p[4]/vol

#variable pxz equal -c_p[5]/vol

#variable pyz equal -c_p[6]/vol

#thermo_style custom step temp etotal pe c_pe press v_press pxx v_pxx pyy v_pyy pzz v_pzz pxy v_pxy pxz v_pxz pyz v_pyz

#thermo_style custom step etotal pe press v_press v_pxx v_pyy v_pzz v_pxy v_pxz v_pyz

#thermo_style custom step temp etotal press v_press



compute 	outsp all property/atom spx spy spz sp fmx fmy fmz

dump		50 all custom 1 dump.equil id type x y z c_outsp[1] c_outsp[2] c_outsp[3] 

#c_outsp[5] c_outsp[6] c_outsp[7]

#dump_modify 1 format line "%d %d %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g" scale yes



#pair_modify compute no



run	 	10000

examples/SPIN/pppm_spin/in.spin.pppm_spin

+6 −4
Original line number Diff line number Diff line
@@ -19,14 +19,15 @@ create_atoms 1 box 


mass		1 58.93



#set 		group all spin/random 31 1.72

set 		group all spin 1.72 0.0 0.0 1.0 

set 		group all spin/random 31 1.72

#set 		group all spin 1.72 0.0 0.0 1.0 

velocity 	all create 100 4928459 rot yes dist gaussian



pair_style 	hybrid/overlay eam/alloy spin/exchange 4.0 spin/dipolar/long 8.0

#pair_style 	hybrid/overlay eam/alloy spin/exchange 4.0 spin/dipolar/long/qsymp 8.0

pair_coeff 	* * eam/alloy ../examples/SPIN/pppm_spin/Co_PurjaPun_2012.eam.alloy Co

pair_coeff 	* * spin/exchange exchange 4.0 0.3593 1.135028015e-05 1.064568567

#pair_coeff 	* * spin/exchange exchange 4.0 0.3593 1.135028015e-05 1.064568567

pair_coeff 	* * spin/exchange exchange 4.0 0.0 1.135028015e-05 1.064568567

#pair_coeff 	* * spin/dipolar/long/qsymp long 8.0

pair_coeff 	* * spin/dipolar/long long 8.0
@@ -34,7 +35,7 @@ neighbor 0.1 bin 
neigh_modify 	every 10 check yes delay 20



kspace_style 	pppm/dipole/spin 1.0e-4

#kspace_modify 	mesh 32 32 32

kspace_modify 	compute yes



#fix 		1 all precession/spin zeeman 1.0 0.0 0.0 1.0

fix 		1 all precession/spin zeeman 0.0 0.0 0.0 1.0
@@ -55,6 +56,7 @@ variable emag equal c_out_mag[5] 
variable 	tmag      equal c_out_mag[6]



thermo_style    custom step time v_magnorm v_emag temp etotal

thermo_modify   format float %20.16g

thermo          10



compute 	outsp all property/atom spx spy spz sp fmx fmy fmz

src/KSPACE/ewald_dipole.cpp

+13 −15
Original line number Diff line number Diff line
@@ -439,8 +439,7 @@ void EwaldDipole::compute(int eflag, int vflag) 


    for (i = 0; i < nlocal; i++) {



      vcik[0] = vcik[1] = vcik[2] = 0.0;

      vcik[3] = vcik[4] = vcik[5] = 0.0;

      for (j = 0; j<6; j++) vcik[j] = 0.0;



      // calculating  re and im of exp(i*k*ri)
@@ -458,29 +457,28 @@ void EwaldDipole::compute(int eflag, int vflag) 


      // compute field for torque calculation



      partial2 = exprl*sfacrl_all[k] + expim*sfacim_all[k];

      tk[i][0] += partial2*eg[k][0];

      tk[i][1] += partial2*eg[k][1];

      tk[i][2] += partial2*eg[k][2];

      partial_peratom = exprl*sfacrl_all[k] + expim*sfacim_all[k];

      tk[i][0] += partial_peratom * eg[k][0];

      tk[i][1] += partial_peratom * eg[k][1];

      tk[i][2] += partial_peratom * eg[k][2];



      // total and per-atom virial correction (dipole only)



      vc[k][0] += vcik[0] = partial2 * mu[i][0] * kx;

      vc[k][1] += vcik[1] = partial2 * mu[i][1] * ky;

      vc[k][2] += vcik[2] = partial2 * mu[i][2] * kz;

      vc[k][3] += vcik[3] = partial2 * mu[i][0] * ky;

      vc[k][4] += vcik[4] = partial2 * mu[i][0] * kz;

      vc[k][5] += vcik[5] = partial2 * mu[i][1] * kz;

      vc[k][0] += vcik[0] = -(partial_peratom * mu[i][0] * eg[k][0]);

      vc[k][1] += vcik[1] = -(partial_peratom * mu[i][1] * eg[k][1]);

      vc[k][2] += vcik[2] = -(partial_peratom * mu[i][2] * eg[k][2]);

      vc[k][3] += vcik[3] = -(partial_peratom * mu[i][0] * eg[k][1]);

      vc[k][4] += vcik[4] = -(partial_peratom * mu[i][0] * eg[k][2]);

      vc[k][5] += vcik[5] = -(partial_peratom * mu[i][1] * eg[k][2]);

      

      // taking re-part of struct_fact x exp(i*k*ri) 

      // (for per-atom energy and virial calc.)



      if (evflag_atom) {

        partial_peratom = exprl*sfacrl_all[k] + expim*sfacim_all[k];

        if (eflag_atom) eatom[i] += muk[k][i]*ug[k]*partial_peratom;

        if (vflag_atom)

          for (j = 0; j < 6; j++)

	    vatom[i][j] += ug[k] * (vg[k][j]*partial_peratom - vcik[j]);

	    vatom[i][j] += (ug[k]*muk[k][i]*vg[k][j]*partial_peratom - vcik[j]);

      }

    }

  }
@@ -517,7 +515,7 @@ void EwaldDipole::compute(int eflag, int vflag) 
    double uk, vk;

    for (k = 0; k < kcount; k++) {

      uk = ug[k] * (sfacrl_all[k]*sfacrl_all[k] + sfacim_all[k]*sfacim_all[k]);

      for (j = 0; j < 6; j++) virial[j] += uk*vg[k][j] + ug[k]*vc[k][j];

      for (j = 0; j < 6; j++) virial[j] += uk*vg[k][j] - vc[k][j];

    }

    for (j = 0; j < 6; j++) virial[j] *= muscale;

  }

CRA Git | Maintained and supported by SUSTech CRA and CCSE