From 8d79db03d38061cf18c12954472ba64a79aa12ad Mon Sep 17 00:00:00 2001 From: julient31 Date: Tue, 21 Aug 2018 13:47:38 -0600 Subject: [PATCH] Commit1 JT 082118 - created pppm_dipole_spin.h/cpp (child-class of pppm_dipole) - improved pair_spin_long.h/cpp - created documentation for pair_spin_long - new 3xN fm_long vector in atom_vec_spin (with associated comm) --- doc/src/Eqs/pair_spin_long_range.jpg | Bin 0 -> 11980 bytes doc/src/Eqs/pair_spin_long_range.tex | 20 + doc/src/Eqs/pair_spin_long_range_force.jpg | Bin 0 -> 16016 bytes doc/src/Eqs/pair_spin_long_range_force.tex | 23 + doc/src/Eqs/pair_spin_long_range_magforce.jpg | Bin 0 -> 9440 bytes doc/src/Eqs/pair_spin_long_range_magforce.tex | 17 + doc/src/pair_spin_long.txt | 84 ++ src/KSPACE/{pppm_spin.cpp => pppm_dipole.cpp} | 1012 ++++++++--------- src/KSPACE/pppm_dipole.h | 213 ++++ src/KSPACE/pppm_dipole_spin.cpp | 750 ++++++++++++ .../{pppm_spin.h => pppm_dipole_spin.h} | 71 +- src/SPIN/atom_vec_spin.cpp | 14 +- src/SPIN/atom_vec_spin.h | 8 +- src/SPIN/pair_spin_exchange.cpp | 3 - src/SPIN/pair_spin_long.cpp | 370 +++--- src/SPIN/pair_spin_long.h | 5 +- 16 files changed, 1839 insertions(+), 751 deletions(-) create mode 100644 doc/src/Eqs/pair_spin_long_range.jpg create mode 100644 doc/src/Eqs/pair_spin_long_range.tex create mode 100644 doc/src/Eqs/pair_spin_long_range_force.jpg create mode 100644 doc/src/Eqs/pair_spin_long_range_force.tex create mode 100644 doc/src/Eqs/pair_spin_long_range_magforce.jpg create mode 100644 doc/src/Eqs/pair_spin_long_range_magforce.tex create mode 100644 doc/src/pair_spin_long.txt rename src/KSPACE/{pppm_spin.cpp => pppm_dipole.cpp} (67%) create mode 100644 src/KSPACE/pppm_dipole.h create mode 100644 src/KSPACE/pppm_dipole_spin.cpp rename src/KSPACE/{pppm_spin.h => pppm_dipole_spin.h} (66%) diff --git a/doc/src/Eqs/pair_spin_long_range.jpg b/doc/src/Eqs/pair_spin_long_range.jpg new file mode 100644 index 0000000000000000000000000000000000000000..bc133d10cf0a0a2a07db95afaf64961911f94084 GIT binary patch literal 11980 zcmb7q1yo$kvhJQ?a35SVxCXc21B1I0Ai*_2fZ#y}2<}dBcXtRH2niZof(J>E5C{@T zkVpP={&UYd@2-2_>*>{dSJhY5)q7U=s_w44#k&oFKv`Z%9sq$rz=!)EaJLM|0w`ec zpL$;)_Z1Zb6%_(O#X?6%!@$A9!NJDD#>T~i65!$y;$dSGkP;9Q5tERR;NX*ylM<6d ziAjk67y+T&>p)O3QBg68aj|iU|G&%KTYwM)cmke70TBXVLJ$fe=xz|820#GHy|?%N z{#{Vf&><*b5XQZf_}|$76yJRUa8N)1m;i+U06@1De{21Da@Z;f{X6}i@92*vfkz+s zn(bjItSs1z|0(`28c*l&Cn{q{F8b0D7Vnu?8WJBWa?0F?{%G4X@r%|WKFapqe}VtJ zfri}k%>H7bnfMVQ`c1jQY%~bA8@BDtz`;7v9b;OtUr|^k)gTjb!M{`~`iqc%f&bM! z2y?UBVT~0X-T1>I@AssO3>l1O9DX7BP`V3=KX65STO$_FFt&&0;7VidI5)E7!>>Ig zO{t~M1-n-$O9yA0NZ-3*$P$-uPXl~Sx~y`#-?&jD01U8rG7Hx0H#g-(BUa0FRTSWtVBH?B%u|4a z4RMh1Lj2|VWb{3rr2JUS<03!l@rtf=Z#y4Ff&ADZR~nB#Hq~#6#&QNRMDn@ZI}y%& z-+%x_x+VZ@u-SxR$aMK8gcWNlCVhhE$p?ediO^qPnnqz#8-57+Iw{RBqq_UBSjhpP zaD{o2X=>rot2V#hETVxCC?p;ld|tA*y! zM-fzpl!NBiCG2bf7%vV}=_55|Q(s5}goBs@%1A(5T?PTz?jmi_a{m~*FA(0aO#9Dv zNC1SUD=jv<2#n$KY>N-Wv_^fV|AF|sMTw{(gZ)ihqkjd!7&$PdB*yDU>cQR*LKUwV zaB=A(PfL6QZMMVCKD_Fc|mK2!sm0MKcCmT#S8X z`F}Y<^z}@!zSl0LHQMseG}$;_is?N=YPM8a_WnYQlvbVkyS^zML6l|42)+%j;7B!n zOt@FS^i9sOP*@CO&(_3u6=zwyWNn_ld5a~!FPI%#iA!A@%0td&;NfLivY;=>kO(fn zg*EO2cIra17Xr2`U9xB&UiU!kxw#jzUVhdO+$~Ed0L?pHaO0l^@$^16@GonfpWtlr zlfDB8E}d(P7}Y8cYg^~m)JI33eQF5y%qhXftIuA2zTZ{#)4I(!+iu(UxJ&94%3QkX ze$(5X6|n~=Qor93U9fZ|;+SYiJ6m?Ur9-IN1D!nJ@t>BY*QGLjTSOyThLP0xAE=0^OqIqZGP7%VAO(&}ml zX7bnPTs7Qa0em)u|6$Mc;UbDDv zd3l-VGZ4W`k|DgaUh-98_i!a8^R`>PF9<*E^0fnz@J~lSb(+*RU*T+UnBe2#L*7@v zlf>U$Cd?kfVh!A}C%%&9Omd6mKZU;?IXDXPSp}bVqJDA59$633yBDx7D$pO*CdXab z1l#}09Y3B3$oLuR-P?(}>uJsuwEB=JMirFvwwLJ9QuVmeS~@Ll>!naG58Y zR{_!E>-;=PPmPfp4pmX3kLKvd zT2U3Sg*nzaAER1U4laH$Zxw56wjkDpfoGI6kbs&NKtzRuX&SklnS@$b+H8OQF_8ev z5ts10LdbS=C(|W^#**NnXv+xnqUfqbVVX6r(I~dqPOz&vEb+yA2&tHJ<%|%*g2vY# zbuPAfNiw|7W87IXG+NChZa>C^%nz60jb|e^j&p+1-~q^Uosqt4HVjt znAu%o)lj$*nZNqwT621~WxOHE#}s-^V7c&ZCi_{0cBi|JM>@&WA{m~`pme8H1W!Ad zon%4&`|=m^kCJ%@Kg!#QM^j&e7niK;tD!~KR$TymGLas&{ezL0s^#8zbvz=LK;Zq; zIoIg%+ss!l+5Ort+eifH)n|;U`3;HAM@d5DuWVoYiA-B0$#CDp#g|0t3STiWmN|LgxRq&fSL+QXq`H{{aF17ef7h z6$RZ#MQG?aXc%ZH_h#-X2qAz9C8FacriW{pyAr_oWP%b4NEmqd1+*=aQX8hG(HI3a z3;U&IbDvCaN0C!k9%L-s^P8TiHdCATLQ`(`+bES1#$c>(6y{FWo9yd1QC z7W(_c@#IcogF#)oEa4ymt?UBZ*M&q{Ps^JwFsR3Pnf~ntuV!;af$i%u^RxH|#G?X{IF`{J^?gff)ahUKkGWaP z`p6z?p>$Q2IdPh%=$iSx|HZ( z(ZB-THsVpcc`*<1ush%z`8azHGA<0Gd@`3Zdtd5ZUgp|(thQ^j6v$;kuZD^m(AxE0~P0m0BdfE3Ob7vtH zuLEQ68spRfTv^0(c2&$>k7&XW9Sxn-p|9@4*g8t?cfjHMksozuje*i0f|a#<)E=y&r%3(^y=-K4kCq9c9X)N$Q$+7hVH!}UcapdJu!O!ON>r^(0axR? zLAl1P1IIVyMgQsUzlJbf>%*j=2$t>141d z@6F8W*0NT^KC9~~A9g!me}3-WDmKJ0_BEkZOV@D<52Z#bV^C1{h%s1iwptvQI*L}9 z)Y=(hA)?#F=|D^>EJ-=GPuWo(yq4=Djom!8I1osMZSmNq#HK1&p+?}_jzjmX*O7_9 zc_iaAsMD9!gnb{X?B8$3>(5?0HieJ%Xh;h`_doP@08Em57Ry4Bi_`R2^BuYaP}TH4GL-5ME*jXB3DiHdAqUa>5C|3xzTb@m-pSd6+s8WyB1|a8Ag61C0`_hFHn}Edc7UI zO0#DgJJ7v>Gq1~q;+nFus5RB*A=!W%%Lgs2cfh3ZtqE^+b#TGMb~#l`*v*R#j-voG zM<-WR*WY5icYvyEE9-{+Tz;L{+(Kv+CHrO3laa_YRp1uKFN^oLXy9;svW=9j#ue`m z$M=VOcL1}bTE+|kLgoFVETZ_f{9e@m3q0u|vjwfJ*Yg);9k+jDUqzbL zY7Z^=Ql$7QJy~Ns6I7{WQz+4ZlV_bhZJ{+ZNe`WMexP#A^8M_xc+2U?!)E!nhTHDk z!PGx};(?cpDN4|O(YWIAr0};A-f9}zeQ;r(3(sK}Hztc@EyFknJbo!pR~D7l$=gAC zM#LVko-;sE7002(c8d~|kAI~rIpzxFk`3uS$tJ0eld<31j&*;?pI6lL^PpFrqJ7bn z!u<1739E}FrQgSHhC7vIqxN$9Mp44{RuU#tF2c)%YO)#za%sh?%p|I7*UQaIhU=a; zmxe-aG>WG7Wf@}J3$ZUBF?~Uufv$WPI@X$WbBM%Mv}{Fm*+0+rBQ|Nxtu=|ZQ2C)P zx6_+5;6g*TH)L&~^rUTcI|0L9y=vxY4fofLVH!<#h%~yLLB{i0*TfO=lCjMPM*9JG z0Gch^kCHJ5%wO}rk?%<1a*VT|RU}l;rqWZmwUL+`yy?`XSikH(C!GtDIBqDWiLY`w z4`w}YJ%rOvD16FtkZXUk@WsPI$z8mYOg|!BhjtU`DEAAGb#_Kq*X{Ky;UGb~UKLMo zzIlS#(KHIGG&dR9Pb*WdeF5ZEm7^Rqd?y&P>BCjp-PhgQ@%nU6vXh*Q)ux=uyohFb z;A?oSp!Fh?fKX4AqoEmV&gE~V#IuVxIT?3Asfb|*iKD^k4ze;{BRD#3W^OdmQ7J%K zv8(O0-5I6uMr5pVkE8X5+qa z!MP6ROgH-;(*ynux@Q&56^0puIZjRk#Amo#k5g_Fo1<$R>$_V_{o9U?S{fD{6@(@f zTcW?FIu3MM3Os;ooZ#s6PQ+?7QRwlN?Q>|dR9N_u1^?1mjd{y%J*Rng)DV5*B-W5N zHJTnpZvT0PLv{1@$0;nY=LhfioQz7R^6xj4Ii8VtWRHx}Gb=3PY&v$8UOsU@b6s_m zRb?-YOEY4Az1gd1P1V`Hbr$oEdA=~rr>M-{H#m&!U@JnK?I09yz@&zU712ieS+vE( z<=Dxgf<9NpSO+IyB4mFvH!V$sCNvU0b*$??-CbJcmI&Lc@3@rcLVh@Om&5YW?${`T z?otFnZCV&|6)@5Z1rKT0kvSEFR+Si@cDsFT@_e0wDRVzUU+55_{!Bnm^*|}}yp^Q9 zBpXS~MbM^-j2OC0lg+)kvg9hKzOu7f=CnQwKl{~??$(|b7o4Q_hPhjV!yR`8J*e@^ zd6;y@|=*%Z#8y@-e7UPBFo& zkU+HDcX=8vFykf)UDsKf5;eV6%VB~_=vposdK6#y2Ip0iXn30ohwjO8o>bts_w8+4!jQQcka!#qp@DR~_XYt}W76+x? zoBjrkyh-w>?(x=~L`enj4b5=liwG&yBu&+8D$+>#x^mLl0opFN0G;oIGEZH1S3ga0 z`v;wGS_S977f2pJ73H4b;xexD`01Uxl#SCeHCtFYAPA3<`#SwR>TtTuf`k36LQ(vc z*7EdAS6P)2?ET=!`$Xf|E(G`*>g*2duM~||{WJ1Q(zmpE>#h}iF7n9uyf1tbsv(HpV};o%Pj7jQe|SL zE}g`Sl?`Jf;t5VmOT&CMs6dyx5*NKkDbL|BZxcxu6cuq@nmm5-8MmYB7usZej??8B zu=cTl*kg5wK7NkZyfSFu=)5d-?8Mr~B$eYjohOSax$3KJ?SB8ZtC5B^ zczf&2qSCe0rz}0|luUl-4p>WaffW`wzODQeg$Z$po;UoQW&2A!JOyDBe}2C9rtL~f z(L%WVmCl%gDD2nSFFzlL{bE8m!R(mp(Y*CzHR{$DCJxo!s{ z%NMnJT)U+Bh`lVzn0PC47XMOToXgjrJc@xqSxY88;NkLTS(y~y+`_z+KF?h1Zy#6h zKSygJbxToounwi1L|+Q{(Lm{nzu}102tBV#+nfRYkd~+Z7xdMxbMM!=%3f~lNmJZ6 zWM8;XPf$&qn%p{>XrOrUipR6F(>h_5{%)vTjuQa~ix$))jfARfdt+Jy<-GIr^UIyW zpe_5-cpDK{=C?oV&p-PtwYCt^zR1fIVLf~>b_uJ{U;0S;O?E=(okr+^P{8d4r?aCX zDcb7dC{sH^1BPZ-)4~uEZNkSn{Ee8>V;;yk=mzcPb``X=39cPUAP)5h3bzu{&Tz zK0%qeqIOCJr{W9$)ml}y=Cu6`=IE-$cf2><(tbd2C!Q98lMyyy?^OJf>y78N>u)Zp z+Vwao>-Xo&CxadWOYW{36rDY#LKl}E=U%VzCc?OXME5ptG!ayaBqQgj=(qDCf4HxE zH7+Kr$eiT^0#rRcS--`3Y>5@E%0oS0$aL9+mwU!B5RMbhQN>E%?j`wf7oek(vDD+} z?s|-tcGj-_d~ZqqL+2`Y=7#x3@{wT`Nt(#j*``;r|9Bvr7O2-RFq}+Ba;NsbFs*Xk z4D&LJcr!{HqpSXm)Xp<(WmSyBncZJ_E0Q^>mIk9UHk&B;RJ~DG*|UbS-6=twY*M0a zkoZKu*bR+oYsNnx#dm#LTsf+~WE=S`)>rV@_$N`IJS<<|D~yrgHRy8G2jS^D>=~lM z9*U!ewcY_P=N5(Lay&_-QdR1mCY|3+pB_4@iTEr!F$%)UlMKE{ylXev(5W>}!}z*n zl{Vqtv$G1+JIN1C zfb&zv8?h4`=Iv|=F(pR$w~wLC_3I0ft13zY-zpzl{m9A>xVYcevlHB({okK`LO_uF z6VE?ee-I%ZzY7(?NW` z--wzWOp(`)tT!v=8*3{`^1iibbwmg)HpB^SIf*VJh$I8D;+>1@Br&fcL+%7DsTxx4 zcYqQsBA!8&dT>oZk#jIo)=uf_Tumi_@a=J{8}=_+eEtuV}tezDDl}qr|QD2W)h1 z)D3($H+022;4Nlgn^GVsI;gxwc5GV%A~59Yx0HY@=-fIszZqKR;*O|Us@La>d4nu+$mExwdb?NFM3cd zyBEp6sBXmv_`FVxGmY4QL3VbL+4@EF`vGSN%LJPD>mb#R5D!RMRs9sg7_3aK-?E-M*OiTNgXoxPAUWQq$X$+7Uu z&Ka>&6!BYg_g=SKeprhyA7akE@143JR`xog5Z#GF`}3_T#0q~W{i+{S>BFxs#n)PP zA&gU6Bioho`044wj#rKMrue#0 zzJ8+Hr@myl10=+qJ};RkpvgG7CO@mHTr^IYDX~S|0s2*m^uLlxGp&Vzym;j!Ip-qD z(DMXlp+Vir4#aCh-M45#2{l#D2=P8ZxN-Y>&}=?l=3R6mPfovx2Iu2X$Es(dh240>(HEj; z=WHi>pQMzT6*0A0`lx;cWV+7()_NJ}v|%uRdk1J_wcLhu3>B1jtE+Dt^23l)5$Cv! zC?BMOcj2tLB122gVI}xRNOoQrQUNE80@`K=}z# zukU~aSWe46w|HD43~y+ksnix%tG~Wkk@QoD90gp?#6aZLh54oAtt4D_r2l7@(bmGL zZL1So3tAFlFt#6IX)ZAzfpeCg_OyYSi^7{aGwlwbE#7x;;==x9emPG0E{F-s?bGGZ z=UN%|8-iM!udj&Z$b&*|XCi!aozVEx1FaYqj|B1Bv7j56#}?w$ny-@~qVzC?wc;dD z7N218u#1>Lf(M&0Q}h6X-bc~NDhRB(cl-@mKtbI9vpgwA&?XUR7=>( z5D?q}mTSs=Q>b%F*Drr`1lQW&K-T{!k<&X#Kz1bEI7y1NJtX^>^R`_Ln%w;oc z*;%jbe@I~})ja(YBDmST;z&AJ*&nKc)X_c>C{`?{qN+1{;`D&azwBVKhHj58>xz^% zWnY3g3N)D9!JLTGFeM|%#Hwzb@Kd6+&lwvVQi^xt6^cc#n5}qICK>g~&%T5zF-R1C z=#M5`Yfi64uM`vKaw+g2fs`Wn(Z$!MS96ji-18cyh7+X4HQ+w7Xk}FQC0qLXBU*_s z#;SvKQskJ&*psA8`S^|kdaxvHCa^S>sX7Jj89s_+@0Rd+exas9xnCIX1k@9RS&=bA z$l6~RW08OI49x*zmv9^UAm5@OpN+_}ZAm>xBo0#%p~`bHVRJHtdU=dpT?&&Jsr;jJ#x|@vgVISj9_-zyFV!(4|_#)nPH8b2+!aUNe8)+ zJn3v7wOA`sqgOKp#MzJpZTyBBE`8iaf#q5u&EjS*IOeh$6^YS_ zpLm(K*T-ef3zfAAiQVVD{Pj2wFOe{Z8Ks$y81p!e&IzY)KUY>#m59+zD?=5Q(LQGd zvU7N6`YPhfW&{&H64Qc(?v{Xi><%zz3>yo%e8rQHMP$|c`a@mml|{*?)T{fcyK}T= zy`mySA5|ThW*W}EW6}*0FJQj|{CwG%eE9Qk`_L)OYngeT?Opcm?S{(S=YSvqsnNju zc<8UZ3E|ZMX#s4>O@Lw+dr%2;}xVrjt zxkd25UT*yZaUUoD36w$iIX|fXNTRvVx}m#|lL=tvJgz~SiT%=M1yjcje_vePTfTvG z*%~K4g{nL;F{AMn{Qfx64taW$y{_3gsr-XgAd6VxQ}p&YjL2F>oQ1;TdAT+-ums(! zM-lFV(l;Cx&FK*RT#?$maG8;cDYx{wPhf)wZ(diB)G5k#jGYDl-f=Rvqb z?@ji7E*r&8k@SwStkn1cqaS}U^!>hgKweXoNj&K=%`{hD%V8MDhb}i0C((c=mTYc_ zmG9Ty?zS~pzlu`BD*j7r17pfiCaBU*Xe_B0nG|J-IHtW(Q;Q(fodhbkglLiMb}~Hp z@AsJ9Uq+yGpuR$t)+CFZ(Z!{wXOWsx3Knc|yY5ZtdnWMVwrfpKR=9xr_p=D{qIwr+ zfm-;BX5y#PBh(|e4=!&`gi~4r0Zu8j?490}z#%BcehErWLf52gLL#1u5k{rE!R7LJ?$k+tP!;01a#n(G5~Ic==xW8Fry} zLih)Am<@$6KfN|o8aoZzrtrf8!HN-(vH>8JnW?BM%D4`X=*s;m$Ma`Ao})AXW-1TS z-l*c@qJ?yB)Nhbjd$NGsJ{YZGaj@Mr)kXgtB)12LW zzIg{Euc4M7%6j+ZP*AkTvUL_G3kY@erbAq}Wm7Kbe&WNYS!~$QADA1+&#K|QGZ1}N zkaC;`!tItRHGgL?WReeL1g_)@^O-n@%mY!8Qd^ z#m0#32-1bw>3PV!%Pfeege90lYXI5`%Ys}U6_ftP3E*KBHo-L5p!TR}lmC32u)%|> zM&qVZ6Poa5u;!}81_y#*@RdMQx(+VDNwr1s_J3Fyt|gyy9gNvloa=goW$2l*G56)H zg*D1gCs&%JCHa+%A7;$sb4%|#ElOXaH1;zVwj?1ZD2W*cd|aAck*R1 zgG-Qp6}$<59t_XF#J3sGA|^CHU~({k;FXppBLprEmzpdN+{-~aX3}=mSSySBwO*VL z7?F7E*6<>F+=7k100ueISxQ}ePL;m42X@2svM?I7g^l5G)Mvxne(P$naJ&tKg878K zi^liV>(A07&FQl*lWkk7E()G2f0? zT5G@3pmX;p^+i`0_@sMOi}a~%-vOtgCKU!**QMBs9Tfi4)i9$>8a4K5O+{Slf}KZi z8XC))p(v{l<>7ORy)FcW+(Fls*)5Is=ch7mFmjZ(O3z-oZq9uSy-7MDI5Hk}! zPg9@1l;+Oh$YvOhw5AF=L_pwZ<1TtkNZ zP0y${KyRg4jSHRh24!#^*5rd4C`SvmA<25_z;r7F`R`#W z-`F5}YA`yl6x=Zc4Jxj%w)ilm2K^Ft9BNuw^p zR~V5Mg3Q)%&Qjp|lr7u#&84?ERT3t8NPBVG>`n>nVBl7wt>_(u8*!7;zm3EXYOyV3 zu}83sh_k*ZSF1$aN5J%g9cMVI3~3ZT z+zP6eg*et2Il18*@5jRwp(hJFPp?SW{G7Alt^=z5P2p9DOQya{WqG`$aw@fShQCwT zjoN}?{IO}FPMtgsT~kjkWr6vl1xBO-H0orN{@PvK;!OcjhhF}btKE8TVwU!ZI<6{O zg9NUyfufPavsuRMV81U(L#Y?%TNhK>$XrbZ1kqVe=O=JJxE@?C_Ktp>P<2%tj3^A$ zjQ0G{0?j4w?coZ8@SBa`k4agssn?`oG1pPJr)!j~v9#qrD^`Z%SibMz7pnZq3qAu0 za%iOB^x^QlX&@+=R{Zk^CiEWxHpmZlym~rehTfEA`2-vb>GqJSk?Re0eYvj>IS+~j z!^sFA^97&%YUx*e5{&QS>c1SPYjT+4wXqeif*dt#cPB zz$-+dWl7c?jGYg|Kw4G9!=D9)x&qQgF>ibt!5uHM7B(&(K0ZDgCLs|49uW>6KHeW8AUGHg z1Rfb49vKfE6&>&2F29EXY$QMvLIelG2Ef=LIBd}GQGg5p0bsa4=Kiz5Bfvor!5}19 zE8f52|3>h86F`H5@nFGW0RZSp;}7-)c3Bn>%R9;Q{hvbq2LKjhi$e)v3fdhsTY;KI zoSX^(s4(H98i%gA$&@T{N@0`&el~fHd;vR=a~(wB z{B(EoYas|pGDd&qat8o#Q^x=xQV;}VJGk){h27z%isWH8FmH+o#UKD+G8VybL#F@$ zT=@XlqtWt=iSw1>!bAYsOZp)41(rnh{=rOEi3^{23kQV8$-5$KMw$^a|lUMponh$k_1TO>>O2G zx^a=;`&5qc7U(A;8&zqjbelXr70Llr9c=F#ZRCQWl=C8`7u zW*{83J>tq&f)T!|?MMM3ErTgMw33u_q+p@#bv_C5uc3;a!$_<0{NM z6bYqV$d2<^ET@>{@rE5@xyk9Ut>`+_TQ^*oBeHE-pJj) ztm(>zJPiOLcmtXIpn1hAE?nSqrNW=K|GI$5A~>l2LBedc2%sRpVx@HQ@0m$4;}@b%GHW+te%`X{zimBW@sW!TM=R6u~*P6vwCF*wB?=GnhW*qDn4x*i=pHkzTh!C3DUs6aS zRqcJI$cy_-%a*Skf!t+EGK?TEB_XFM&4-NCgeDp{t8ngI_sl->Vt)R@(lLWn^V=cR zAv~QEL*~tR^~njXSw)a>j6gNl?nD)22Y*{5;C)-q8a^x42Wd+pCp6r~4k3*?cF`#M zOtK}(4x8JN4{RngB-T;s;V9xd)34K->>;h0zFKvR=r1;qS9^Lz(&}Pev75}gxK!Ro zXoQ!q-n~_n`qt?kiyk7z<8u38q&)ogE3pDWrsf7I^<-ze;JmzT<&{y0W>~p&OQdfO zedS$CBX``sfDVBoJLemm@89o?ZnIZc9ylbJhc!y{rTQlS#S~pKMX2MdNW7x(sD8WN zV*$gIYuma{SjsW#G=#|6urhs`i?WI$0=BV7*&Nb3~}FosYK^O`WdMEtNXnz0xGA~t&#pms(>}pqph(1WVy;0 z+IcavQcG{pCrsYg?aeJRWJxEuhwVv9PreZ{q>b0Xes7{Zakm*J3#nQuWSN(-CO;nI zI!Q>mp7$hiQ-E^l9$2aka`sC-i=9zxD>K6<2(Z78zT+`9BdrstLiAVLU+v`&TD{s# z@xK&@Kaxisg^*~8x?EPUe|vkRZ7xgJIG7yvE{eA*YrP|hS1sY~RiVL_wT}W+&1`CF zM|Lfpe_<}*^2N9k8m(s}K5$+`=93w*NK#;^q1w*`^nL*uHEFIbqFZi+nrJF9-d+v= zdUujr*=c%ld5P&7Y2;Ce^+Io3=S6u(nqxOt5Iza}p@5b=s^o|P&YlnU?Gk3fb-XpN zz#F`ou^*Jax5oZ6C)O z=&&J|8oDo8xoEh9Mru1V->*yi#upO8B_9&Wkvcc6E%7TRr@S(^Um2@8IkbjwKT7a_ z*I%r|=E?_Ye^T9J)%1&T_?z_Zs4k!G^a^*%q`~t2Ie3g|e1}a9QJkd`9Sj`uu437Q zvVewR#m=d9bY>tDnxK}>^oLwJ#?fYKRc3c1g=v5x(3kR;=1-TlcFRbDj#4pn69!z0uJSssa{;10k|k+FnRXtJ?>w! z4Hk0VttQxbZ$ca!9cF*Vt`sF`Rb6WD!?id=f0&qXUWy7znv zorkWps=fL*5SLm#$%xI()}{jThW2YqhE?ljSIbE@^;{&U;c%OGq;eG>m(^6JYuoe~ zUGzhWmX`CVg#!%$L_*q{iPigOiSKbQq4B|$-C5JM`(*4E-D>I$th?X){KN=S?Yu{c@5kM(+O&9cT(#Z8Hn=1$xF6^ z9;Re?P{F?eYWHu+p?1ojw-nibxt*{lZL#q`3rr9EGYiK4m)Y@O<6F^RgBWbuEC%8g z`!6WK|D;gE0>B_J7#;}@6&W4@5heu$hJ#=OIPkdC+!9!LG(75NuEA82#R#-q&-gUl zlGB=7@aZ&7%|lX#mZUsWONJSErQMe=_{Fs>tadN|K8PZSfxk;an~2_*;D;ynHQ&BJ z^5K-G}Pdrp+8qGRxQl2QmsOul#wn+#Y`( zXTpszUO$Ts$^)78k9lTZpU|8fq|wL-Hx9Y$>+$relvIH$N_a-G7PmRTRE#nM6mvP3z42D5iit z`g5aWa_3oAwoaKD$BhUYkMk)P^ zHe^<^ijeMNSe3|(^vr1QJ$VRJS{95@1Aev;jy*SuDCyjCftVK}zC zi?R;Q=HqKJmkB6Wm0x16aFC&NErh%{Z7tBRVdj&0ZOf|JOEp`#JcXZ(GtYi``DOJ__wU0j~uUKTS#8sO||Y`H(p^q-EDX(h+C9Ii>KfuFKrpiR=D?fYXsN4tbWw}nr$$g&r2e|TyLSk2 zsVq`jS5DTGlP<(^wy|CAdW1h^ne5T)p*c+YV)*#!X0KwA z;urk@c+d$B{EP(|8^Mj>yi&A{gP7 zb{Th#boDll+E9-MHvBZpIg+YW-lP|#*|FBy?no05&WW(6}t=@*k1`;3jXzedVBNH`B4yWWtpT!Wswfmxsj4aY*X{! zuTV&(T9}0m^BCsUD{cJ0zKyeX>p^AO-{Y3k)kh8~k^TaA67as@m#(FECLlb&z6hZk& zvUX6t%qWHFLBcg&jViaYYS`H77=Lk#m#RpoIe(R&F8MVPziA`dc#|M%pf^kTQS|A) zpD`n+FkdCsBjRq}mm88gGWXd7?b74{)9O-JRe~v@-)aVHA+6Jf{H~vuBB6c>*Eg>! z#PqVw^~Rbu{48()K06)JJ{L`H-qUYyTU1x)>HWl8ncF~qXOYN9G#msA$Iwnnn$ADo zgH=hIy@w}eL8;3A;mD+_PP@@XkC%!va>bg`a(Gogzv)P^96mojWH`9blTG5Ku{K*ud+cxTISZMc4sD-{ zpgj`mcXcYLtqnNYttc+}S*(Tx#(MDYIvr2o)v7Hwt4NyCF1>aP%y!)vz`wOu<&UjP z$Tv@we776OJwdEq%*FC!XP_Tp_izCJRA?qSGRm`q++Vi<-WlC)Pd~SoB=dz{ z?TU|-%uIRKRBffv($^v~akTAy$L2D|Y?^X>jT3?t+{G}1Px>ESS$iax2=7yo2jm9& z?+5iZ34$fw#omm#{PIg*T#&T=vddasD0;iyBVd~1fJkQRlkzQcabKpPO>_3()VJ5tY9p8d$hKQcwmIi@nB z}JwGt9GEav=sG-_9eEQw1T$wR1L=Ss(q)QVG zGq3iZa7GJe+0+$DwIRr)TDgc$Kdrq`5U#w?tP8cE7MYqo?2j*IRoVH(gz^c#pr7U4 znIkhs=FT1y@AX|#%>^htI_+PviI)&Cl_NZP1L*a)1Ud-Vr`NqD~n`0xkPN^MmZK|#f3Txs?0K5 zO^!O6%xZsHTUrXgPHv1UEum+rIaos^U3KjYz=)5OK&r7)YSDc!V%2&X5N)n%O~6zh zUlRg}uTc)|40t9C5BEkG;hAa;UGq%!cE%H>_5bI4iC$xqQ5;@^OUmp-L|mU6-~27v z|K(in0`;#@y!IrUwUx~m1vH8BmHZ8cw=rj1_AY1RnGF^1CFoW_?07MF;^Mo`KQ^M@ zs2$){x6TngUpNN8%4D-@wpG@zTP zAubCb`?Y9smFbG=jph<3e*1JXDoZs(V94Qrw;R_`oY}O$NAnAm7?t=Pv64?)$})m} zb^@=xCW5bkn@V%|Hvsn=7}OL?w+wxP+w|v8ZgA|Dy%Pb(vyU(fq*PVp`^Ub|g)SC7 zGG3C|u$(PT)0KLi)AIOg?RBJxrV5#eRq!{qp-BnAOHZ|XSMQ0)q9Y_WR3>CKJ6fy| zWICU3{xZ??vC^jt!RHwuNV|>y;T+z*z+1Mev|n2wl~;Aqb-1+lg=Iopa;y}m*n1%g z;$@!G1Jx>*TK=St$89OY=u4Uc9do!_P+M}hk=R)#8(6>^@r_ENUvL>uUXJwidno-4 zq?NKz3cp-fvfC6@>06YeO7jfbZgUz_Nvlm#eJK?6tz?D-Q`=QLFlJe;w)!{lA<>L- z#NtOK9$)gjt^IMK$TO>15u`otlv*gAi_eSTpN8B%#my|2|XZLiNUIAf-2?W`Nk`BKnlthUlMTh9mKqjT(66=|F> z8Rdg_UE}G@h^a~{W0J-9ls`P}XH+C;%#dAXO&RZfH^wFR7vZJmdsAE?v|JgmeDJb= zJtSV2oelspH_Wy|^}O6S4r@@?i93AR&p0Jhcuth3_K?$5aLDb875@tpkeGf1rLnAM+G9YSexP zO>&AR{J5lQ*0Ps^)TzEarB3VmD!GpHC;N;Ry;1qsiX?xWg9NCVQ(cVK!(L8-YGN6? zGSw;5XrPjs?r*?C9lIpy=evH(-+9zN z5{BgLQ8%Gz`s7BNy)T+P_mHP0lD1dw`3c~GfP{EXY@5meGx$qYenhhKVp zk*`Dx42t@Kb6>T5NWU2EE=xg_FpwmtoPE7V7?_q9r|&0LF<_!3B%)zr0RNvrZP{7MtS7FBj`;RY3` z400^!uOadKpbhtm|4MVJ%%gbwTA{(8yhiSgSt&>tW&ABaX9-*3P!BPgNuyM4OmEzWjXx2-|)V_4=W50 zAeHa@$;ZnR{KL&kIPmrP#dzfBCij@gN*NhpJco}HL8CuBWJ;$sq!z#M<}((Ad~g8- za4r;CP|hX6^h}>ln;Yz!6N8EZp#*m_L7@3oLPY%LuxfJ!YL?vv8hq`hbh+(W`HT~7 zyWg)F#je0+ippYfmdI~qJ6hbGtySJSL+AFxBM+p}&c~K|wO5kMPmcm&iyL2rHH9gS zhSQWNJ>&4u;oDj7o)8ba8>OCL8+?DWLZhv}PAX-quE!{g>fJW2fxdodx=P1?@Q@=GoOP=p$=$ z>)wy%)eXxm2ZrBFw!E(A`M-ddXqYgYF#wQm#9Vc1^IBI%@jQU$W@!w z7v9Q$JtLy|u2oG*laF5lmOt&mq*TZ#>1eMlD@U)P2Rk~R`YLz)?|uv0l&yA~ovDcs z{-jI(>&mN;WBzX=-h)sizZe}X(M=KcOpwUf3?-6yzh_L7DtV}9u_?zSu}Ic5|t?;seoBfo)Xr8^n< zkB)R{%P(WlE4DL)MXB*zdx*Df4^p_G39t>}rf)M4eCl z?bN=X-A%c5)%ZuGl$e8e>nyVSg8lj5Fji@Kw{zk^EjO@$kESJ zR~BaEkCJ?%r5XypVXn!3*~0NPaIwN&H@Q%O_GSH&+C`j6A~v>1wBzv|y>v)M4{LBO zG7}qK)1ASQqQvUbpn7vg2d6$cTv1Y_y>OyFw7eOgVTpg@a&O-di8PuWdM)dT!`bvm zHI4ROqxdi+k~<(?~Q$4ca!;3jMF= z9LUxp7T%wShc&?KBYsD!Hi83f`gBQrRspQ+_Fr1sPQ0Q!UOY6#`-V0`LHSaagPALX ze(RVrIVI?MKz=op$>nU(tqJ8df>|PFV|< zsXHDwhNES8sVFCZ+P}fEF-r8&Uf&BQ%f^*3q~(yJ*A%;cVDkiTCPP=j$tt0GjlH4p zBP8>xX)`!#BJyq5;zgy;t(#|33=Z9UPLir zA%pUlR(~(+7SFzLBDoS=tr;XBtJ)y8Oove+k050REF?aE_q?268f8!?D>58>;I#-8 zqKW*Opv_`Nznk7p6r9QMMBzFbyutm(KfE_{vBYxazbP0;moY6+r{ z#fUSWiCikCa`|sY#(HXUQsl#e>~3&LK!6U`2e`|tF@C*K6ZvBn$!}?- zGhLmUX<+Gh&G{+9h$A2K=a5(M?+iWn&N zODvW6%-v5_&;q|@wA@#P&f&Zof;;#^OOGf|ajY^yQDcdPA$@vV_bqF;%KTs#TfJmE zW0AAw-sbYF38T-;jH&JKx#Qd`@fFbrt}A%`lK%I)r(!uzk=HMqe%aqq&IOD-G_D8o z{lNS@VD|Im^coQd2TFcmH$u~I&TIiSRq;Jo@ZpxWPieNVC16#& zUa6VA_V&{pPQdwF7TWpawo{pRj7H{##B`H7%gAvU&s0gpTgk00+p^pHYXG zhh=^6riY`^+E#xm(shE*CYRmy+=maAZwC_s{N>N(!-yH!HyojoDY3cgv_lXqR2E{h zUj5oWfu1a6LrLK}#lXQ0ItvKjOMJ1j#Ktb>nQKB=v|Qk$c2Lj4Sy2O*WFowB8Mm^L zz_64`gO7zMTY-?7viE}>w4TrLbwm*$z4bnRz3Gl)Jml0}x$J)#bL>t64x z&bsbXUuEaYy7MloYjX-#*Jq1f@BhPc=vuY?oGP3pPowE|^n{qLCKo zHYIUTdCHF1<`k#G`4aLfdVQHWL?AoU)3%0WSCCTA%rQMRu{2w9D=3yMGQPkWmKFX zIi$(LF#USu#l7eEF3)0>W#R+8un#XZiksAKnEX>z4zYW$?#e@haIIo&Ay#tAY0?{U zn*Ea%udwzo3IiHf0XRZk8tj674eywVlA$^y^om`8^n;fceIy^wwt>twk}vLG z#jtu*NB1x&u$p=sTlm|nNY*mGF!Vb$7^7(8Eg7i>G+wc0gAp0t&^htCerbKNVlni& zRJrv-|7s~=Cf|QZM4jZf1_x}*OKXV`OA%9(S=NbkrMmJAO*hmF-x|0+i28MJS7gM; zEBEGK$}5i?vJixHzV`l%*09=kDk~|;=cvkSe-E~us$-+yFt)R?ZRYRnM(@9^s zdAip?qPDi{abaP0n+5hF-0v)b4#jOxJr;u^1lowz4AHSF7P^^aO<%_Pgyw-%S{8y2 zKK^C)SW9N&0qc*1Iogr0I)%5d&egDYL&&%e_A0)mXq=$-YZc4U_JjvLkHm8+GH-^W z%gmqK`DX+!?p>7DN$ViK9X6FMzm%={{BhDFecFT>r}bg6nAsJ}l+}{e+1dPM31Ub`;*qo3K3rtT*$TLl!C)d@CEz$>fgKKWBlb7E9h>?AH*+JQYS3h{qa20CnymO zYd4Y+|G7s>rk?RtIf@|p2R2ZzKxsG}evG^w$7^5B)7D>+gyG~bs{jO!25DN+@v=+& zM43>%E>h)Vec4peq;JA~#iki!fB)5$h1Oca&EW~ypTAw6yD4I0<3>aH59d^$kC^Cw z#1zA4LV^lD`X)Ll70e`xDkI8o8~b60QQV=n+2xZNS#Qa!I;Q`MzzZLkTIA542K@jr zI7;G|NDivwJu;DVhFE=r$9!Jo0rg~Kgp%Iq5fna|PDWY9N*s^=QZZ{$bX0s66KX@p z?44ksiy0dxV)WKG2%V2y0?U6K+JA!+=Xs7RJN1!~r4PWua`4ky+`u0fZVJxQF4k+E zr90i-3xbV#h_j|heF!74HSSWfN)yxAoUJl?vXMm(5^?y;ByjUqK0SVxR+uOCw>Ph{ULaGcOtSJLJu15D+{w@Lm6=_&y&16$M z#P$eCdqzhn^JmOu(_PC^om!wlx#-udo?XFnS;_=`1kgTw45}`6JWz~I=HQ*G>fd>D z9Ezi)kt+bDR%a`lG=|jqnP$hH>m%8>96MX~k{hw561w>>OW77RvDg%eQev^#7{0;T zvs8@dvE%Hb@Sz$_#jBXdk=IgW<{^||rNlLIfSRK=NrU)%7whyt<9NATiW4m&!Br7q zOZnjF9ikX2IXoBBR|;n4Pl>oeD+jy@^7gMmV$^zXY!#!%9FDW_qJmcJh{ z2H;@FB(PGTbSV;af9ikZqQdz@{I~kQWBel=l7a-Jh*KrW10gU*dDsPksrw`3F8~Jr zLBeV%T?iBcApgMwe+Z}m@jnDE z_%9YPG9(D55g7nOsc^)pu*B(-B2Ea7_`AUTa zQvn5VkRUKj5{L?R68n$YFc?Y)My3M)!_=U^kT84-GE6m0>|Zp0;J>@{j}`z(Li$VB zKZ*bEz+YT{HNkp~^M9xLf5?G30r~&Y^l!=eKZpzei?cU>PUm3#x2JRdM!|#yTN#7F zFz18e{(9IoZkcW4x?kn?g~LlhEd-i1%QVYDjIBlyRLc0uT#vYRGa{1VV< zYaWP)0=Bn#!doK2WF#jaJ`#Z@fuJq_Lml(Y1cCiQB5uTv+5;IAdP5|p$Ct@<<~f#3 zY!1{B@iT)K%|(T9jXJbMg5CPYaDKl53dsnD0WFzTsdJ*^$Y{^fI#tRrCSotdC;IW$ zMJqSu;uM!k$q)iO8ZjT{^|jirTCO;x?h;R=k9No({1OW*9t7~vdfmF>2 zG_I6b{My`iLN{j88fJR}XEeB9Qiw zlggG-s#FRJ3W4bbGpU3iIOyJla@c}wUwyE*F>&xAhRM9@7|w2IIads)*fVYZ1L>&p`na?2AEVi`cSNDf~%h6;ZGm~&n6DvV(*Gkp18 zajr&ksZ#f9Pb40vS(-QcZ2&FN5lI64p`Y)o2^BgT=Hi$rX@6+RPgmMxi0+iqw;WtA zX~zwe-W^NN6b)a{V2%7Iy9e$zQy~YBvhx>hG~al80AUPThcG4yW<*b#Z2P->KWT5J`eP5$HaaxH($|MrmTNP1I8z32aV2Vfm1Lr#(MQFPQ8YbOOEkXbr z5q=jJdR8R;?w;HBhbUWZLTnSZ3l$X&)uCk_L^}vvgalbSqkjvDXpcbUTIgRL zVEM&?!^F&ALP8h4YU-C!0v5TT^<=*!D!1o0CIjIKZ-eBir9BD~MT93eJflO9WLOay z+J(qs9uF!}4RIB?MBfoPpA3u91Wagt+wE{Blf$U1Q#D=9fU1svK}XTF#|4Av zje;W?>SreLN+v8U6};Hv{+uA*h#1MY98a~G9MFUG&*C#EWgh1=>iZS!#w}RaU(<0F z3DKlJqnQ^X1l#W(ULKQ1%*B8fd8;h3V(O=HLtBfxs_Fw4A67R>z6!FBAU&?Lbx=oK_hs_=31rU5SIyN`vN3q`C19HB2 zKB3jrusAe%w5-EUD@9|Xenfe`wTX3)R|l(w}LC*y7$qMTeg*zL&_k2t{GghCRlS89DF}EurB(gwPHk64I z%&Kgm;|L-Z=Gl%LFG9|8Au$X?8Sf2GQkk5y+^L&VEX6Ad9A~;vezOv^HN6K7OZ*h2 zhXue-D|uE#DTBlW;d!ltDNlr&BkvrF|t(K96_((awn`d*N!O>fU(jvwO*}c>9 zG&cRbaM@|Qk6W3su|cwVffj+4L6@&tsmLNg2Qm?G9g|!VsM2_K6yviqdlqWCL!~6& zD?AkZ#=a`=%~lZ%hg8za8Xi3~LMyOxGzIuuH%%E{nZlf(H+Y&}$Gr;OpR-QkQeS3?c^bO`56PPxYY7&L!G-(Vk-4?CwGZH@*nLC| z930{Ui1~~MH!)LXB?>48YXw0x3MlakgD|&;E)m-u>H|ACIRJ^x73+Qww5USLA}y-O z?_>jyYH@Cny=icJsW~jGYtvW-f@!qpaPLWQfOU(Y`>Ti#nl^P z3J@rqbeRmNc`IZfrMtbS-P!}a>+WX|NlAvLnLPT)KjkhvMc_#&PM@7n`BBC=PEt#w zn^x)0h`Uao#mkEu2+{loY>X}NO{D>h<)4U>&j+tV=A*=!?T6-1(A^r2qZ8q1hvT@N zid5n>2!!=~>`HCHHZ{B?@b0dtrxvC8k-aQ(!QEJoBvqvQ&WpD3Z z6gA$lm9cp>g=&=9rGogiJi=77bj*-gvL$cuv1P)Zn-A@T0hn?jzX5l;+319srW*#E z=X1XSoTKT-pP{FT!36;mNq*F+dv082oigftTw^o!Q36^-JG`nxUR-VRRsxM;4NztR zUE491$e=Pi8+gfZX^i85Ks^aWai6!j)E0B%3QcYHv7veWRnPFE#PRfVfcV0V5O*@L=jMR>?>j#k3cLq1{2D5112vN4EvfsQJ zU(_KMFVaT;yi!yOQLs;4Sy;cq4CC=rkwxP_u3@Mk>h^3#Ls1+?t&Ag-ZffgatyHDO3E8xs-H^g+w*(r*6Xm+T4 zjk@Rjc2RrLAVUs^Zx80?gJM!LF#=;C;&9^RzTehg+}&FLj0y&oE|g4?NIUfIX9_mv zQNIojD%b@CL7nH|@R{1mAO-6+?bSht7d?z$_sdIHQJNeHmWygksOVh_J0dWeMo#TY z`iD371rFJP{c}~^*QcXZd4Y(ajilbF9h@$SReEe3u1T)PR1Dde$}`pFQNBU! zh!}}EkcPk(W;7VojkMk?*GUEuOfHO;`uXeHY|TZ`InhweOIys8#HwG=_q;_1qq}za(tusRdU? zVmQm6nNxp~mP~^O1^9J@(HTlBgLi_nEI-Uf%R>Rq(i&^?dHQ3wfr)fhskLgfO0{t* zUqe!}fd{%?94#TCRJBuQdrY@#OB7!9Xb%i*yqD!j;-H!tg*?C48mRtg0R#CeBNHwu zZE7M_ucORgIpiD3<>mB6c8fFAC|)Zgel<7-v-a$Xdv0fMmjSqfjJym00)c>D^aI>30FnR%4E~eQ z21O?r77PZ3!f-G#F|hD(@bGYPaB=Yoh~fAIL@ETh@L`k(S@p@NJ~o@qXS1~vqdOFmwths}eMrwxsQVjl&c z7;z82H%#n_2HT>tUJPm5NkmDdDjh*5t?$qv=5hK*{qTyrFUIEZ3jpF1QvQHxa%w>KP$aC<1% zPKvru7_U69hs~q`@4i#s2ktk4aIFKLQLr$oorc6^z=p((t`#t%8eI^h%9`&&(2z5x zV@q(&8Wq8%gSwc8v5AY!lDPosSk$KGP=Ixe&7W4_9LVrso75TLCIFL|pIM}d0`@Sf z9!6wefivVtS}8F+7@kV#0-CsvA6`QI(hUIs-s6PqMT((|f@9S~K_>Kg5E}r?v=gi} zh*AGw-E&<9&{TZ&NBVCDldAp81xPEl0GQ@4iDc+Mp?`5Oqf5iX03d*8aQR=*f4a1U z5SD*_{M|thA|xC=!oZ+EqaO_U8w!O%FhF1cPDYPOL`qD;Mb03mN^y@7!H*tXSm+@K zf`P9Rhxyu%UUhn76Hb#PN6GCQ2P^K}0t>(@Idh2vMBejNylv@>p;aCMgNgO#K5x9| zTBMMCQ9^A;&B+?=mCMJq8^ss=94iKOrJor}5$<^6U}JN8+agm=5P_i>iZ6pro)5yXF=rfMu}74%6B zFU6HKyI@1o5sdiAR;6X=D38wY9m)#fST*G2jj!#xocOn!+Kbp_@Sc23F4CHzHnqK> z`My9o%-^$K>yIoX|3s~Ia|^UE#mKuJ)O)2CvKKBF1y%<7WJwqY(INM6@5UM6yx#uJ z&Ok{Cggb4vRfCXU4~_*LmJe7?VJ_}c>Yt8LSepxO2g%NAF%-FPk&4rt*O!=Nsf<3GCAM=oe<&btR~B@Uh5Pu5hM`yG+Wj zaKW!s&JDcZU95(;CQ|!>H_>OoxX_l6xyD7m?Xg*4%W+5uL~?6E1gK+S^*dGnVP z*?J*11dPfC#p0A!Zs?z5Ida5LlsnSa!=bs{(sxmt?3&KtRbp#ym_rO}ZS`F|Awr7a zHG3li*&?^JS|u;HFaG4)vok+ERHWFQlX2I6+)}@#?wt&-S;15N8`C^ z9%&DVya{W?*-m%^(jR==zdT-i_TF9AiQ&}9vuC;%6q;>hp=_(D}0?{@p52o~& z8zy%qo~Iu3w|o3@3!b3V&F7fAwzGC~^CfxqjAffr?fqTJPBu?6WUv4oUbMoDgTpM> z8U{f|Mm+4CjzLN;-%P%>vT{yA|1ebzapjNn_j;Gu#uAb~nnOHw7V~;TOs_)v36IG^ zBF~}u^^?+1xl>2S7tJ{E6{~UX9G-LD3qqp)aC;a(Vf}a8^X?77XR{R@<3dCU1d13k=bJE-+Dl zn*$&qkSK_+<$qy-|H)85&o3Y_7zV?KU}3}lOfcxx3rYkK!$=so#YxqS9epw2^o&eA zd=l!>vE^-C4?iPRO`K3M8DnI;{9=+Cn&T@{rhXOs|CoI+M8TI5@@14W8D##EF{VFD ztCxF}=~;{Elz(%tHhDQx+zZb3yrC-1q~aQb8E*G0(JEy6Ps zQl+mS5jB3|MLnn_3&wEhv)~eBvM}ObSo)Tmp)UI%KjLH1DZ$JQ{1quKs0-Es?5dD)E6`^ z)~y~4DMfQj=wbZXww8G;s25jNE2##`REQ}H zB6nXpR=q7nX5X}kc}iiXeRDbu4bzyGJkv;za$Ehn)7v(as(i)vdu!FIz+73q0Kb{L z|DA^~V@?pS0x7?WvCYY9Mw|dKgJ{&;r8#RXJna>L}r6H z`+~h?XKT^H@$q%;R&_S|o43IA(0fv@Y^mdYUos*7ZX^Db#q}j{5<$A)_s?nmmO%* zQ$)CuR&iw{rt)?dt44oaSq_&c=UD;9JNbi)l=)s zeB14ceEt9#OnLp#*3ga)F8h&ac!Q(rY;(-}_^_wnDzLIWDnOdEX?$!5_`7j9v95Fi^y|U>8U0VW%NgjmH5`m{WY!I%KLb?;WwzBmP@Z_lFiCFHIn#SWgX{l4C!zcdA6zX$27XTrHv zUpG5{eL=;p-k>QTJy+oL#ZaEzZu&=u)9awzffmh`;uZOfe8ZfG5AzM5MVxOGKBW<} zM9EU&6%^JH;&?D6=BG)q-BbQC5l({tU|qjD<^>aT;YZsev6GtJt$pFd;uTuAtYaHk zbL#ed6J6DLA@$=7BM0wk8X=C5?Y`^{zLU~4EAR#+QPbH(C?U3w1Mi16+`%@t`)*}2 zA?~@pe!^~LNOaH^tIJ2h?1)wVq}W7KXWLhVjbvK*yi-%$&yB`Y5A{1NmeXe!GU2N1 z!X3m|NLc;n}BO!pb8p@-^mqYV9L_z#h}EjG z{hwX-En(52^$y4VwU9DtfB~2dZ(+GEC2B|U4i{5!v z)K8j>P|r>S`=vD2Z5x%y^%yAxu1dsR%errl#2Oz4GV?jjVl2bgA1`{Tg>z(X_HG9C zSg9KDyEX-Vr^9>AEHs&!|7p)!<(mq6GhXCj!lq!yL@Y^JHv8?U@|_4f#ZK*E8@uil z{!SOkTU+jXa)B|H{h;;H^ART7`fZB2iKXuq|5J-}*_FsP&;lKF^LN_eLUhNvG%oyB&Ab8?QU7Y1K)7CW-kK_#IP2 zgC1+iRCRkcn}}Vc%ApsIe#rcE2`ZL9FBZGHw){K)+)DfL)R%Z>+kVqwx?5o2Uz8jl zm*Y)6SWg#k%CZaLmX&C9<(^G*TYY#1~&dM_dy(E+pf= z?v8;+UG=5oSA4bLHItl(kA}KYKEFGP6E4d$`O0~DLR_?E$?1J~K4&X9a5tI(R#8sO6|$Lej=BBg7aTdqQlPAD=txChG%9NXRhMR9IAR`sdw#YUve z_`4PFI2?^O8XV>RPG8w~*74!fX%S)4)LgT}o9hv6wlAYAN_LP#+8-@v9EE%e;{)qV z+(U}ZgSB*2DXeMFTI6bd5tBZAv<4Pw)C-7Dn-5+PEnM$wYwJoXJ1WBVYQq9^lap&N zr$j`gh5Mk@^1v&*jmIm~u$tdm>gt^qg;GK(jU7aB*YN@&p{mln3XE#Z?b`y;6XD!h z-jmZj4BzFh#1BQ&p6gK}EFV7PU*^t;BPX61mzc60sCpKldkf(8bOfkqsHu!u)MW3z zIa+12cA4s+>@PfDF9iDC;2z6G@GLgBdA{1sQuH~34UUd0{4ECYSF{Tq21B0@)Ewc) zzVy*;`&{L!pT))wR{n}Xh%#ISh!L0IT+=J~jJoU+G{-jc{65=}#1{|aLhw~*;=Ct3 z^9WHn*wKp38Z{TX1s+i6es&A3O0dY$uv?e14;A9_AHyoG6XgZML4c-xn*TiKxN%aV za&0tG1F~2mi*l0a&5%6ByEq5WSDMI-7G_&KipE=jZE#Oge3qC=tG~r}$JSVSj0I?j zVs*k21}Fu+Iqce`Xb62}+V$v++Ti{ZbKDl@eE zY#VgL#I;Zb*8cEe;8C=oIP5W2|Kum|;!k=lD!f5pET~O(4xeQ2m1pL^vE-w>d=2VM zQ_W_U1wiQ^#&!c%XZ08lAo|j6yi@c!%~|(L1clyA8+`cQA)6b$s^&L7=_Ik@dlOix zY8S+6Ml7$IC?qS#N3KeYIn051(Nyl|2L7G>oioO~DIsXpkF@Z+*L2&r6xZV%F??Ib zocw&E56__~Z5y=`{ZHGxwjPgVjfO3G=n4}`_F>7gBzWC#S%qqYx5qxfqxr7N9t)73 z+NIb%HS4Pfah-9i_IdRv=69AEJo)TBaf0$75t1g0?j^g+@-RtW*N#i}Lci&a#UL&} zk40x~PO@OxHx{i6?XEq5d)BfjLM-oGJY`yM97l+p`y!)I?@RqJoumPd#KvE+8h&Xu zzEJX%MyDQks~!V&O3~fVrOwivZ)M{(1e@ehkC{f~>|g8+oK5StcFuoFzsu}EPb!d8 z@Z=>>QqpI89zF7(oJ+QO{JD)XWk#%S(wk@wPsiI8Cw|2eHPh@ia);_uuhHcCC(k*GO=h&ERaLoD!2G7TKLdjP#Uc_|F7>ban{G=j5Tro7ReB zvL6yk6_-+Kv?7{XDH7jaUz1`vFiN+>Y46aXWV|VuH}3an@}6>6Hkd2)!Qw*I_ho zKT!6Zc<8c}cV8CD6mKrPaXqP=71Cz0NE>gZp8T^NbCSt2-PTe^1;|8q^Qld%8QK@r(f>OI}@<>ETf*#E&H16eEY1zhXgw_lh~QcglWo2Zwb;(6v{9T12O8ec`ZsF zr@YMAd$My~dXEiiQuatMk+;C-t`%uSBB64}@=UnA0?IC>!T}XFDz{WpttPb2 zSY)?eXGJX*y`$|dn79Qf=KXu5fR}?d^aDEaG(szW9Lcqjjk+ zv&rRH2|uD7Hp>VkKY#0}3Lp`o)q*nSH=Q=u%1y;4v~1QShfIL zuX>bjb7N%DJBsslwrosj?nd$L%0mGV>GQ6#j%_;wKnosVG(|SdL&}dLr%4=CPp!`~ zVLH6=PIU8BfGU)XJyrF(Zcux@%HYR@$0cPK$Q^>v}E6# z483Ij?M?{z6W)iT1-`G$;ae!DJ z@WTesf_~WODxf~te{_iDplBrOAC3GYLJa&Vhz9>7{UQG!Ae!TY7W|)Zbo2gQLG*u# zV*ZJ9!vC*v&c8F5;pnRo^i2x-66Mc>w!iOS;A-^7T#mkCs?lTn2mjo?z^{T5BU||k z_KkZ$N)~RR#k;GPUKH|M%o(p0DA*ELJvP=Kah-ByfR)aOkE|>f0jIU=FRufZ!-Li7 z8Mlw#^TZS;6p6Fgx-zac>K?#T1kMNQfFR%%rXJnrIwlj%AgqQj%Y#|lo9u7}YFfDl zgcsu(1rc63=c9^Zt#{1iPU2s6D)0J024jpYSrxR? zNBSagtcf5l8M2{=d|?_&{;=>5Tx0bfcy7ETTdw;a5;+5$!>UPnvgZlr0Y3>k+xw?f z-_dcfO$1$oSBjRvKxHgLwHP~>I3ZXt-u7|6%an2M1Ab5A<8&jN^wgf-^DBqS6n(8+ zq6=7uUNaQ*l10+mN830^=3d{PeVUHf6n1GIv1N6KxC1N7b5YMpa9J7SV04r@Aad{S z;{ASm+`zAQnWnlJi8H~5IIn4huYG0fjw3o{n7O6L4Esj>ovsoNG;o@F;Bs+Td>*!*Ze(3e znU}$ot<1OVQ;#fy>@lD8n=;%2ZTnPSnB$TK046@smdFjgk~tQrzYcRS6;e(&(x2TF zV-=h1V9nNGstDaD9L>NP8ncfWA~C*?o-8fVf_dsmp9Bp!N zsAk5jn};MN%SRBINc3M)sQPPyocfkRUak6B7UXZOWd=tg40(Z@1gwbG;hT$v3Y z+mq^6F(y&9gV&?^*H<)OQv#3Gyc3{)q_XnMBzM3;hXP|CzG{5>l0=w^HTp))d7(K$ zXu!UwQ}g0Zo8u?&@J8cZaa~jrJsz3Rn*Suq1=}P80kDW5}u>#}k=e%+KyvKa- z^s&8FG;at-N;#)r(LN$0>Nl=VZA2C@M@a5Z-1PLu+t0=6_50jxN*m4?ets@V^1w1p z1a=t2PC)Nd>nBsHIQ=cFM=H5Gn6zeh_I}h^0)7>sd?)fO_(h<=vyVr@t0>CR8Puq&zS{3R+Gx5dvRgzp~I=ZR%lOu;E zGR^rM>;}3zC5MD~<4CvZ9z34M-^3l0cDOsuh#>^96}Fl;&n%u-N?>d;g_HRuFtVEA z`1lR@OyB>Uv;6~CS=qz>$;>VAMv}%6Ap;PV%Dum%x%z3wEJav07qU|!5lW>IAJwR9 zpe>^Px*%$6UY~Us3>-4`ZDfj}(%}%FC-y(X3aUgCLH@j$^bw|M26FT1=Ky~=(>)b>$xw`3-xU^ z*NswszwapT-X*n_-ty8{t8xgZ$9$eW8V~zfHvFq(mk-$}X!ubf^9~n(o_h$BfK$3X zFC<)$2q#jFE&WMI#W0t3F`!8#v!WPOlqmlCt^RTrXV?2XN#heGK56`Nel?8LK4Y?S zoeNQ|Yn*L%5*T$pCp@te{m{qwPb2=-2>bGT4oDem>e3* zvCqKj`O@{c5BYp?7EO}^(t`z>3s94hjgNAZ86^cQ`Tt(?>W4kU-sj4F zb=H(zSCwg9f>#o$7?aEP@Z|1EVo?yC9%d#BLOl}t#e+?t%#>XN(zI1G6#nWNM$keg zXx|BkBU+4lqAw#SHcuq%N!$}Sab})M*R$Vu#2>h)U-)NZfqYsSh*%PSSe(e=g!z$; zl8}6TA>Rvf2eg2h?&o2>8>(7BdcDOA=)r7;%}GjT=O)v9f6FT$%!uNm3Fm1C{xlgT z@UJ_wSDjT;D8>Q=yu>FVZTmt+47U!Phf)6A>$1M7LbiIEc`YQC`mn6%^&7ZJ?4lnj zgtOw!a!?ts?J=pjkG~ZZdFKND%6U`GM^aH-K#567Gsee+VNz-$^dSb17*0wZAKMiv z8CEBtr46br_3;Js3Lle`(wE^560#xoKdQ@C6ZfPb)3G!_FbRMyoq*&OrHn6_&BlhT zHD+aQE*aTjG`BCs4=x(r!Usv?wQQtLP2$XDM_g;vp~D5bA*>REkL1)SuiRzG_JTAm z+3=99A6=s$7*p1SEYgRNn8?gqVBWi4Z@U#+bG literal 0 HcmV?d00001 diff --git a/doc/src/Eqs/pair_spin_long_range_magforce.tex b/doc/src/Eqs/pair_spin_long_range_magforce.tex new file mode 100644 index 0000000000..ad40cf9d8b --- /dev/null +++ b/doc/src/Eqs/pair_spin_long_range_magforce.tex @@ -0,0 +1,17 @@ +\documentclass[preview]{standalone} +\usepackage{varwidth} +\usepackage[utf8x]{inputenc} +\usepackage{amsmath,amssymb,graphics,bm,setspace} + +\begin{document} +\begin{varwidth}{50in} + \begin{equation} + \bm{\omega}_i = + \frac{\mu_0 (\mu_B)^2}{4\pi\hbar}\sum_{j} + \frac{g_i g_j}{r_{ij}^3} + \, \Big( + 3\,(\bm{e}_{ij}\cdot\bm{s}_{j})\bm{e}_{ij} + -\bm{s}_{j} \Big) \nonumber + \end{equation} +\end{varwidth} +\end{document} diff --git a/doc/src/pair_spin_long.txt b/doc/src/pair_spin_long.txt new file mode 100644 index 0000000000..c5b4a7b33e --- /dev/null +++ b/doc/src/pair_spin_long.txt @@ -0,0 +1,84 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +pair_style spin/long command :h3 + +[Syntax:] + +pair_style spin/long cutoff (cutoff) + +cutoff = global cutoff pair (distance in metal units) :ulb,l +:ule + +[Examples:] + +pair_style spin/long 10.0 +pair_coeff * * long 10.0 +pair_coeff 2 3 long 8.0 :pre + +[Description:] + +Style {pair/spin/long} computes interactions between pairs of particles +that each have a magnetic spin. + +:c,image(Eqs/pair_spin_long_range.jpg) + +where si and sj are two magnetic spins of two particles with Lande factors +gi and gj respectively, eij = (ri - rj)/|ri-rj| is the unit vector between +sites i and j, mu0 the vacuum permeability, muB the Bohr magneton (muB = +5.788 eV/T in metal units). + +Style {pair/spin/long} computes magnetic precession vectors: + +:c,image(Eqs/pair_spin_long_range_magforce.jpg) + +with h the Planck constant (in metal units), and a mechanical force: + +:c,image(Eqs/pair_spin_long_range_force.jpg) + + +The following coefficient must be defined for each pair of atoms +types via the "pair_coeff"_pair_coeff.html command as in the examples +above, or in the data file or restart files read by the +"read_data"_read_data.html or "read_restart"_read_restart.html +commands, or by mixing as described below: + +rc (distance units) :ul + +with rc is the radius cutoff of the short-range component of the +long-range interaction (see "(Cerda)"_#Cerda1 for more +explanation). + +:line + +[Restrictions:] + +The {pair/spin/long} style is part of the SPIN package. It is only +enabled if LAMMPS was built with that package. See the +"Making LAMMPS"_Section_start.html#start_3 section for more info. + +The {pair/spin/long} style computes the short-range component of +the dipole-dipole interaction. The functions evaluating the +long-range component are part of the KSPACE package. +They can be enabled only if LAMMPS was built with that package. + +[Related commands:] + +"atom_style spin"_atom_style.html, "pair_coeff"_pair_coeff.html, + +[Default:] none + +:line + +:link(Tranchida6) +[(Tranchida)] Tranchida, Plimpton, Thibaudeau and Thompson, +Journal of Computational Physics, (2018). + +:link(Cerda1) +[(Cerda)] Cerda, Ballenegger, Lenz, and Holm, J Chem Phys, 129(23), +234104 (2008). diff --git a/src/KSPACE/pppm_spin.cpp b/src/KSPACE/pppm_dipole.cpp similarity index 67% rename from src/KSPACE/pppm_spin.cpp rename to src/KSPACE/pppm_dipole.cpp index 9b59f9cd7b..a03f5b9980 100644 --- a/src/KSPACE/pppm_spin.cpp +++ b/src/KSPACE/pppm_dipole.cpp @@ -21,7 +21,7 @@ #include #include #include -#include "pppm_spin.h" +#include "pppm_dipole.h" #include "atom.h" #include "comm.h" #include "gridcomm.h" @@ -50,8 +50,8 @@ using namespace MathSpecial; #define SMALL 0.00001 #define EPS_HOC 1.0e-7 -enum{REVERSE_SP}; -enum{FORWARD_SP,FORWARD_SP_PERATOM}; +enum{REVERSE_MU}; +enum{FORWARD_MU,FORWARD_MU_PERATOM}; #ifdef FFT_SINGLE #define ZEROF 0.0f @@ -63,34 +63,34 @@ enum{FORWARD_SP,FORWARD_SP_PERATOM}; /* ---------------------------------------------------------------------- */ -PPPMSpin::PPPMSpin(LAMMPS *lmp, int narg, char **arg) : PPPM(lmp, narg, arg), - densityx_brick_spin(NULL), densityy_brick_spin(NULL), - densityz_brick_spin(NULL), ux_brick_spin(NULL), - uy_brick_spin(NULL), uz_brick_spin(NULL), vdxx_brick_spin(NULL), - vdxy_brick_spin(NULL), vdyy_brick_spin(NULL), - vdxz_brick_spin(NULL), vdyz_brick_spin(NULL), - vdzz_brick_spin(NULL), v0x_brick_spin(NULL), v1x_brick_spin(NULL), - v2x_brick_spin(NULL), v3x_brick_spin(NULL), v4x_brick_spin(NULL), - v5x_brick_spin(NULL), v0y_brick_spin(NULL), v1y_brick_spin(NULL), - v2y_brick_spin(NULL), v3y_brick_spin(NULL), v4y_brick_spin(NULL), - v5y_brick_spin(NULL), v0z_brick_spin(NULL), v1z_brick_spin(NULL), - v2z_brick_spin(NULL), v3z_brick_spin(NULL), v4z_brick_spin(NULL), - v5z_brick_spin(NULL), work3(NULL), work4(NULL), - densityx_fft_spin(NULL), densityy_fft_spin(NULL), - densityz_fft_spin(NULL) +PPPMDipole::PPPMDipole(LAMMPS *lmp, int narg, char **arg) : PPPM(lmp, narg, arg), + densityx_brick_dipole(NULL), densityy_brick_dipole(NULL), + densityz_brick_dipole(NULL), ux_brick_dipole(NULL), + uy_brick_dipole(NULL), uz_brick_dipole(NULL), vdxx_brick_dipole(NULL), + vdxy_brick_dipole(NULL), vdyy_brick_dipole(NULL), + vdxz_brick_dipole(NULL), vdyz_brick_dipole(NULL), + vdzz_brick_dipole(NULL), v0x_brick_dipole(NULL), v1x_brick_dipole(NULL), + v2x_brick_dipole(NULL), v3x_brick_dipole(NULL), v4x_brick_dipole(NULL), + v5x_brick_dipole(NULL), v0y_brick_dipole(NULL), v1y_brick_dipole(NULL), + v2y_brick_dipole(NULL), v3y_brick_dipole(NULL), v4y_brick_dipole(NULL), + v5y_brick_dipole(NULL), v0z_brick_dipole(NULL), v1z_brick_dipole(NULL), + v2z_brick_dipole(NULL), v3z_brick_dipole(NULL), v4z_brick_dipole(NULL), + v5z_brick_dipole(NULL), work3(NULL), work4(NULL), + densityx_fft_dipole(NULL), densityy_fft_dipole(NULL), + densityz_fft_dipole(NULL) { - spinflag = 1; + dipoleflag = 1; group_group_enable = 0; - cg_spin = NULL; - cg_peratom_spin = NULL; + cg_dipole = NULL; + cg_peratom_dipole = NULL; } /* ---------------------------------------------------------------------- free all memory ------------------------------------------------------------------------- */ -PPPMSpin::~PPPMSpin() +PPPMDipole::~PPPMDipole() { if (copymode) return; @@ -99,23 +99,26 @@ PPPMSpin::~PPPMSpin() fft1 = NULL; fft2 = NULL; remap = NULL; - cg_spin = NULL; + cg_dipole = NULL; } /* ---------------------------------------------------------------------- called once before run ------------------------------------------------------------------------- */ -void PPPMSpin::init() +void PPPMDipole::init() { if (me == 0) { - if (screen) fprintf(screen,"PPPMSpin initialization ...\n"); - if (logfile) fprintf(logfile,"PPPMSpin initialization ...\n"); + if (screen) fprintf(screen,"PPPMDipole initialization ...\n"); + if (logfile) fprintf(logfile,"PPPMDipole initialization ...\n"); } // error check - spinflag = atom->sp?1:0; + dipoleflag = atom->mu?1:0; + qsum_qsq(0); + if (dipoleflag && q2) + error->all(FLERR,"Cannot (yet) uses charges with Kspace style PPPMDipole"); triclinic_check(); @@ -123,30 +126,30 @@ void PPPMSpin::init() error->all(FLERR,"Must redefine kspace_style after changing to triclinic box"); if (domain->dimension == 2) error->all(FLERR, - "Cannot use PPPMSpin with 2d simulation"); + "Cannot use PPPMDipole with 2d simulation"); if (comm->style != 0) - error->universe_all(FLERR,"PPPMSpin can only currently be used with " + error->universe_all(FLERR,"PPPMDipole can only currently be used with " "comm_style brick"); - if (!atom->sp) error->all(FLERR,"Kspace style requires atom attribute sp"); + if (!atom->mu) error->all(FLERR,"Kspace style requires atom attribute mu"); - if (atom->sp && differentiation_flag == 1) error->all(FLERR,"Cannot (yet) use kspace_modify diff" - " ad with spins"); + if (atom->mu && differentiation_flag == 1) error->all(FLERR,"Cannot (yet) use kspace_modify diff" + " ad with dipoles"); - if (spinflag && strcmp(update->unit_style,"metal") != 0) - error->all(FLERR,"'metal' units have to be used with spins"); + if (dipoleflag && strcmp(update->unit_style,"electron") == 0) + error->all(FLERR,"Cannot (yet) use 'electron' units with dipoles"); if (slabflag == 0 && domain->nonperiodic > 0) - error->all(FLERR,"Cannot use nonperiodic boundaries with PPPMSpin"); + error->all(FLERR,"Cannot use nonperiodic boundaries with PPPMDipole"); if (slabflag) { if (domain->xperiodic != 1 || domain->yperiodic != 1 || domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1) - error->all(FLERR,"Incorrect boundaries with slab PPPMSpin"); + error->all(FLERR,"Incorrect boundaries with slab PPPMDipole"); } if (order < 2 || order > MAXORDER) { char str[128]; - sprintf(str,"PPPMSpin order cannot be < 2 or > than %d",MAXORDER); + sprintf(str,"PPPMDipole order cannot be < 2 or > than %d",MAXORDER); error->all(FLERR,str); } @@ -154,7 +157,7 @@ void PPPMSpin::init() triclinic = domain->triclinic; if (triclinic) - error->all(FLERR,"Cannot yet use triclinic cells with PPPMSpin"); + error->all(FLERR,"Cannot yet use triclinic cells with PPPMDipole"); pair_check(); @@ -167,21 +170,17 @@ void PPPMSpin::init() // kspace TIP4P not yet supported if (tip4pflag) - error->all(FLERR,"Cannot yet use TIP4P with PPPMSpin"); + error->all(FLERR,"Cannot yet use TIP4P with PPPMDipole"); + + // compute qsum & qsqsum and warn if not charge-neutral scale = 1.0; - hbar = force->hplanck/MY_2PI; // eV/(rad.THz) - mub = 5.78901e-5; // in eV/T - mu_0 = 1.2566370614e-6; // in T.m/A - mub2mu0 = mub * mub * mu_0; // in eV - mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz - spsum_spsq(); + qqrd2e = force->qqrd2e; + musum_musq(); natoms_original = atom->natoms; // set accuracy (force units) from accuracy_relative or accuracy_absolute - // is two_charge_force still relevant for spin systems? - if (accuracy_absolute >= 0.0) accuracy = accuracy_absolute; else accuracy = accuracy_relative * two_charge_force; @@ -203,11 +202,11 @@ void PPPMSpin::init() while (order >= minorder) { if (iteration && me == 0) - error->warning(FLERR,"Reducing PPPMSpin order b/c stencil extends " + error->warning(FLERR,"Reducing PPPMDipole order b/c stencil extends " "beyond nearest neighbor processor"); compute_gf_denom(); - set_grid_global(); + set_grid_global(mu2); set_grid_local(); if (overlap_allowed) break; @@ -224,9 +223,9 @@ void PPPMSpin::init() iteration++; } - if (order < minorder) error->all(FLERR,"PPPMSpin order < minimum allowed order"); + if (order < minorder) error->all(FLERR,"PPPMDipole order < minimum allowed order"); if (!overlap_allowed && cgtmp->ghost_overlap()) - error->all(FLERR,"PPPMSpin grid stencil extends " + error->all(FLERR,"PPPMDipole grid stencil extends " "beyond nearest neighbor processor"); if (cgtmp) delete cgtmp; @@ -236,7 +235,7 @@ void PPPMSpin::init() // calculate the final accuracy - double estimated_accuracy = final_accuracy_spin(); + double estimated_accuracy = final_accuracy_dipole(mu2); // print stats @@ -282,10 +281,10 @@ void PPPMSpin::init() // don't invoke allocate peratom(), will be allocated when needed allocate(); - cg_spin->ghost_notify(); - cg_spin->setup(); + cg_dipole->ghost_notify(); + cg_dipole->setup(); - // pre-compute Green's function denominator expansion + // pre-compute Green's function denomiator expansion // pre-compute 1d charge distribution coefficients compute_gf_denom(); @@ -293,27 +292,27 @@ void PPPMSpin::init() } /* ---------------------------------------------------------------------- - adjust PPPMSpin coeffs, called initially and whenever volume has changed + adjust PPPMDipole coeffs, called initially and whenever volume has changed ------------------------------------------------------------------------- */ -void PPPMSpin::setup() +void PPPMDipole::setup() { // perform some checks to avoid illegal boundaries with read_data if (slabflag == 0 && domain->nonperiodic > 0) - error->all(FLERR,"Cannot use nonperiodic boundaries with PPPMSpin"); + error->all(FLERR,"Cannot use nonperiodic boundaries with PPPMDipole"); if (slabflag) { if (domain->xperiodic != 1 || domain->yperiodic != 1 || domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1) - error->all(FLERR,"Incorrect boundaries with slab PPPMSpin"); + error->all(FLERR,"Incorrect boundaries with slab PPPMDipole"); } int i,j,k,n; double *prd; // volume-dependent factors - // adjust z dimension for 2d slab PPPMSpin - // z dimension for 3d PPPMSpin is zprd since slab_volfactor = 1.0 + // adjust z dimension for 2d slab PPPMDipole + // z dimension for 3d PPPMDipole is zprd since slab_volfactor = 1.0 prd = domain->prd; double xprd = prd[0]; @@ -381,7 +380,7 @@ void PPPMSpin::setup() } } - compute_gf_spin(); + compute_gf_dipole(); } /* ---------------------------------------------------------------------- @@ -389,7 +388,7 @@ void PPPMSpin::setup() called by fix balance b/c it changed sizes of processor sub-domains ------------------------------------------------------------------------- */ -void PPPMSpin::setup_grid() +void PPPMDipole::setup_grid() { // free all arrays previously allocated @@ -406,11 +405,11 @@ void PPPMSpin::setup_grid() allocate(); - cg_spin->ghost_notify(); - if (overlap_allowed == 0 && cg_spin->ghost_overlap()) - error->all(FLERR,"PPPMSpin grid stencil extends " + cg_dipole->ghost_notify(); + if (overlap_allowed == 0 && cg_dipole->ghost_overlap()) + error->all(FLERR,"PPPMDipole grid stencil extends " "beyond nearest neighbor processor"); - cg_spin->setup(); + cg_dipole->setup(); // pre-compute Green's function denomiator expansion // pre-compute 1d charge distribution coefficients @@ -424,10 +423,10 @@ void PPPMSpin::setup_grid() } /* ---------------------------------------------------------------------- - compute the PPPMSpin long-range force, energy, virial + compute the PPPMDipole long-range force, energy, virial ------------------------------------------------------------------------- */ -void PPPMSpin::compute(int eflag, int vflag) +void PPPMDipole::compute(int eflag, int vflag) { int i,j; @@ -440,20 +439,20 @@ void PPPMSpin::compute(int eflag, int vflag) if (evflag_atom && !peratom_allocate_flag) { allocate_peratom(); - cg_peratom_spin->ghost_notify(); - cg_peratom_spin->setup(); + cg_peratom_dipole->ghost_notify(); + cg_peratom_dipole->setup(); } // if atom count has changed, update qsum and qsqsum if (atom->natoms != natoms_original) { - spsum_spsq(); + musum_musq(); natoms_original = atom->natoms; } - // return if there are no spins + // return if there are no dipoles - if (spsqsum == 0.0) return; + if (musqsum == 0.0) return; // convert atoms from box to lamda coords @@ -464,51 +463,51 @@ void PPPMSpin::compute(int eflag, int vflag) if (atom->nmax > nmax) { memory->destroy(part2grid); nmax = atom->nmax; - memory->create(part2grid,nmax,3,"pppm_spin:part2grid"); + memory->create(part2grid,nmax,3,"pppm_dipole:part2grid"); } // find grid points for all my particles - // map my particle charge onto my local 3d on-grid density + // map my particle charge onto my local 3d density grid particle_map(); - make_rho_spin(); + make_rho_dipole(); // all procs communicate density values from their ghost cells // to fully sum contribution in their 3d bricks // remap from 3d decomposition to FFT decomposition - cg_spin->reverse_comm(this,REVERSE_SP); - brick2fft_spin(); + cg_dipole->reverse_comm(this,REVERSE_MU); + brick2fft_dipole(); // compute potential gradient on my FFT grid and // portion of e_long on this proc's FFT grid // return gradients (electric fields) in 3d brick decomposition // also performs per-atom calculations via poisson_peratom() - poisson_ik_spin(); + poisson_ik_dipole(); // all procs communicate E-field values // to fill ghost cells surrounding their 3d bricks - cg_spin->forward_comm(this,FORWARD_SP); + cg_dipole->forward_comm(this,FORWARD_MU); // extra per-atom energy/virial communication if (evflag_atom) { - cg_peratom_spin->forward_comm(this,FORWARD_SP_PERATOM); + cg_peratom_dipole->forward_comm(this,FORWARD_MU_PERATOM); } // calculate the force on my particles - fieldforce_ik_spin(); + fieldforce_ik_dipole(); // extra per-atom energy/virial communication - if (evflag_atom) fieldforce_peratom_spin(); + if (evflag_atom) fieldforce_peratom_dipole(); // sum global energy across procs and add in volume-dependent term - const double spscale = mub2mu0 * scale; + const double qscale = qqrd2e * scale; const double g3 = g_ewald*g_ewald*g_ewald; if (eflag_global) { @@ -517,8 +516,8 @@ void PPPMSpin::compute(int eflag, int vflag) energy = energy_all; energy *= 0.5*volume; - energy -= spsqsum*2.0*g3/3.0/MY_PIS; - energy *= spscale; + energy -= musqsum*2.0*g3/3.0/MY_PIS; + energy *= qscale; } // sum global virial across procs @@ -526,32 +525,29 @@ void PPPMSpin::compute(int eflag, int vflag) if (vflag_global) { double virial_all[6]; MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < 6; i++) virial[i] = 0.5*spscale*volume*virial_all[i]; + for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i]; } // per-atom energy/virial // energy includes self-energy correction if (evflag_atom) { - double **sp = atom->sp; - double spx,spy,spz; + double *q = atom->q; + double **mu = atom->mu; int nlocal = atom->nlocal; int ntotal = nlocal; if (eflag_atom) { for (i = 0; i < nlocal; i++) { - spx = sp[i][0]*sp[i][3]; - spy = sp[i][1]*sp[i][3]; - spz = sp[i][2]*sp[i][3]; eatom[i] *= 0.5; - eatom[i] -= (spx*spx + spy*spy + spz*spz)*2.0*g3/3.0/MY_PIS; - eatom[i] *= spscale; + eatom[i] -= (mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2])*2.0*g3/3.0/MY_PIS; + eatom[i] *= qscale; } } if (vflag_atom) { for (i = 0; i < ntotal; i++) - for (j = 0; j < 6; j++) vatom[i][j] *= 0.5*spscale; + for (j = 0; j < 6; j++) vatom[i][j] *= 0.5*qscale; } } @@ -564,59 +560,59 @@ void PPPMSpin::compute(int eflag, int vflag) allocate memory that depends on # of K-vectors and order ------------------------------------------------------------------------- */ -void PPPMSpin::allocate() +void PPPMDipole::allocate() { - memory->create3d_offset(densityx_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:densityx_brick_spin"); - memory->create3d_offset(densityy_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:densityy_brick_spin"); - memory->create3d_offset(densityz_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:densityz_brick_spin"); + memory->create3d_offset(densityx_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:densityx_brick_dipole"); + memory->create3d_offset(densityy_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:densityy_brick_dipole"); + memory->create3d_offset(densityz_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:densityz_brick_dipole"); - memory->create(densityx_fft_spin,nfft_both,"pppm_spin:densityy_fft_spin"); - memory->create(densityy_fft_spin,nfft_both,"pppm_spin:densityy_fft_spin"); - memory->create(densityz_fft_spin,nfft_both,"pppm_spin:densityz_fft_spin"); + memory->create(densityx_fft_dipole,nfft_both,"pppm_dipole:densityy_fft_dipole"); + memory->create(densityy_fft_dipole,nfft_both,"pppm_dipole:densityy_fft_dipole"); + memory->create(densityz_fft_dipole,nfft_both,"pppm_dipole:densityz_fft_dipole"); - memory->create(greensfn,nfft_both,"pppm_spin:greensfn"); - memory->create(work1,2*nfft_both,"pppm_spin:work1"); - memory->create(work2,2*nfft_both,"pppm_spin:work2"); - memory->create(work3,2*nfft_both,"pppm_spin:work3"); - memory->create(work4,2*nfft_both,"pppm_spin:work4"); - memory->create(vg,nfft_both,6,"pppm_spin:vg"); + memory->create(greensfn,nfft_both,"pppm_dipole:greensfn"); + memory->create(work1,2*nfft_both,"pppm_dipole:work1"); + memory->create(work2,2*nfft_both,"pppm_dipole:work2"); + memory->create(work3,2*nfft_both,"pppm_dipole:work3"); + memory->create(work4,2*nfft_both,"pppm_dipole:work4"); + memory->create(vg,nfft_both,6,"pppm_dipole:vg"); - memory->create1d_offset(fkx,nxlo_fft,nxhi_fft,"pppm_spin:fkx"); - memory->create1d_offset(fky,nylo_fft,nyhi_fft,"pppm_spin:fky"); - memory->create1d_offset(fkz,nzlo_fft,nzhi_fft,"pppm_spin:fkz"); + memory->create1d_offset(fkx,nxlo_fft,nxhi_fft,"pppm_dipole:fkx"); + memory->create1d_offset(fky,nylo_fft,nyhi_fft,"pppm_dipole:fky"); + memory->create1d_offset(fkz,nzlo_fft,nzhi_fft,"pppm_dipole:fkz"); - memory->create3d_offset(ux_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:ux_brick_spin"); - memory->create3d_offset(uy_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:uy_brick_spin"); - memory->create3d_offset(uz_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:uz_brick_spin"); + memory->create3d_offset(ux_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:ux_brick_dipole"); + memory->create3d_offset(uy_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:uy_brick_dipole"); + memory->create3d_offset(uz_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:uz_brick_dipole"); - memory->create3d_offset(vdxx_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:vdxx_brick_spin"); - memory->create3d_offset(vdxy_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:vdxy_brick_spin"); - memory->create3d_offset(vdyy_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:vdyy_brick_spin"); - memory->create3d_offset(vdxz_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:vdxz_brick_spin"); - memory->create3d_offset(vdyz_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:vdyz_brick_spin"); - memory->create3d_offset(vdzz_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:vdzz_brick_spin"); + memory->create3d_offset(vdxx_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:vdxx_brick_dipole"); + memory->create3d_offset(vdxy_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:vdxy_brick_dipole"); + memory->create3d_offset(vdyy_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:vdyy_brick_dipole"); + memory->create3d_offset(vdxz_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:vdxz_brick_dipole"); + memory->create3d_offset(vdyz_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:vdyz_brick_dipole"); + memory->create3d_offset(vdzz_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:vdzz_brick_dipole"); // summation coeffs order_allocated = order; - memory->create(gf_b,order,"pppm_spin:gf_b"); - memory->create2d_offset(rho1d,3,-order/2,order/2,"pppm_spin:rho1d"); - memory->create2d_offset(drho1d,3,-order/2,order/2,"pppm_spin:drho1d"); - memory->create2d_offset(rho_coeff,order,(1-order)/2,order/2,"pppm_spin:rho_coeff"); + memory->create(gf_b,order,"pppm_dipole:gf_b"); + memory->create2d_offset(rho1d,3,-order/2,order/2,"pppm_dipole:rho1d"); + memory->create2d_offset(drho1d,3,-order/2,order/2,"pppm_dipole:drho1d"); + memory->create2d_offset(rho_coeff,order,(1-order)/2,order/2,"pppm_dipole:rho_coeff"); memory->create2d_offset(drho_coeff,order,(1-order)/2,order/2, - "pppm_spin:drho_coeff"); + "pppm_dipole:drho_coeff"); // create 2 FFTs and a Remap // 1st FFT keeps data in FFT decompostion @@ -644,7 +640,7 @@ void PPPMSpin::allocate() int (*procneigh)[2] = comm->procneigh; - cg_spin = new GridComm(lmp,world,9,3, + cg_dipole = new GridComm(lmp,world,9,3, nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out, procneigh[0][0],procneigh[0][1],procneigh[1][0], @@ -655,26 +651,26 @@ void PPPMSpin::allocate() deallocate memory that depends on # of K-vectors and order ------------------------------------------------------------------------- */ -void PPPMSpin::deallocate() +void PPPMDipole::deallocate() { - memory->destroy3d_offset(densityx_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(densityy_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(densityz_brick_spin,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(densityx_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(densityy_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(densityz_brick_dipole,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(ux_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(uy_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(uz_brick_spin,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(ux_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(uy_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(uz_brick_dipole,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(vdxx_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(vdxy_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(vdyy_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(vdxz_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(vdyz_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(vdzz_brick_spin,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(vdxx_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(vdxy_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(vdyy_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(vdxz_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(vdyz_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(vdzz_brick_dipole,nzlo_out,nylo_out,nxlo_out); - memory->destroy(densityx_fft_spin); - memory->destroy(densityy_fft_spin); - memory->destroy(densityz_fft_spin); + memory->destroy(densityx_fft_dipole); + memory->destroy(densityy_fft_dipole); + memory->destroy(densityz_fft_dipole); memory->destroy(greensfn); memory->destroy(work1); @@ -696,61 +692,61 @@ void PPPMSpin::deallocate() delete fft1; delete fft2; delete remap; - delete cg_spin; + delete cg_dipole; } /* ---------------------------------------------------------------------- allocate per-atom memory that depends on # of K-vectors and order ------------------------------------------------------------------------- */ -void PPPMSpin::allocate_peratom() +void PPPMDipole::allocate_peratom() { peratom_allocate_flag = 1; - memory->create3d_offset(v0x_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v0x_brick_spin"); - memory->create3d_offset(v1x_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v1x_brick_spin"); - memory->create3d_offset(v2x_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v2x_brick_spin"); - memory->create3d_offset(v3x_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v3x_brick_spin"); - memory->create3d_offset(v4x_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v4x_brick_spin"); - memory->create3d_offset(v5x_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v5x_brick_spin"); + memory->create3d_offset(v0x_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v0x_brick_dipole"); + memory->create3d_offset(v1x_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v1x_brick_dipole"); + memory->create3d_offset(v2x_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v2x_brick_dipole"); + memory->create3d_offset(v3x_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v3x_brick_dipole"); + memory->create3d_offset(v4x_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v4x_brick_dipole"); + memory->create3d_offset(v5x_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v5x_brick_dipole"); - memory->create3d_offset(v0y_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v0y_brick_spin"); - memory->create3d_offset(v1y_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v1y_brick_spin"); - memory->create3d_offset(v2y_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v2y_brick_spin"); - memory->create3d_offset(v3y_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v3y_brick_spin"); - memory->create3d_offset(v4y_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v4y_brick_spin"); - memory->create3d_offset(v5y_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v5y_brick_spin"); + memory->create3d_offset(v0y_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v0y_brick_dipole"); + memory->create3d_offset(v1y_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v1y_brick_dipole"); + memory->create3d_offset(v2y_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v2y_brick_dipole"); + memory->create3d_offset(v3y_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v3y_brick_dipole"); + memory->create3d_offset(v4y_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v4y_brick_dipole"); + memory->create3d_offset(v5y_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v5y_brick_dipole"); - memory->create3d_offset(v0z_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v0z_brick_spin"); - memory->create3d_offset(v1z_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v1z_brick_spin"); - memory->create3d_offset(v2z_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v2z_brick_spin"); - memory->create3d_offset(v3z_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v3z_brick_spin"); - memory->create3d_offset(v4z_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v4z_brick_spin"); - memory->create3d_offset(v5z_brick_spin,nzlo_out,nzhi_out,nylo_out,nyhi_out, - nxlo_out,nxhi_out,"pppm_spin:v5z_brick_spin"); + memory->create3d_offset(v0z_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v0z_brick_dipole"); + memory->create3d_offset(v1z_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v1z_brick_dipole"); + memory->create3d_offset(v2z_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v2z_brick_dipole"); + memory->create3d_offset(v3z_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v3z_brick_dipole"); + memory->create3d_offset(v4z_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v4z_brick_dipole"); + memory->create3d_offset(v5z_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out, + nxlo_out,nxhi_out,"pppm_dipole:v5z_brick_dipole"); // create ghost grid object for rho and electric field communication int (*procneigh)[2] = comm->procneigh; - cg_peratom_spin = + cg_peratom_dipole = new GridComm(lmp,world,18,1, nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out, @@ -762,44 +758,44 @@ void PPPMSpin::allocate_peratom() deallocate per-atom memory that depends on # of K-vectors and order ------------------------------------------------------------------------- */ -void PPPMSpin::deallocate_peratom() +void PPPMDipole::deallocate_peratom() { peratom_allocate_flag = 0; - memory->destroy3d_offset(v0x_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v1x_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v2x_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v3x_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v4x_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v5x_brick_spin,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v0x_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v1x_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v2x_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v3x_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v4x_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v5x_brick_dipole,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v0y_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v1y_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v2y_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v3y_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v4y_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v5y_brick_spin,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v0y_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v1y_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v2y_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v3y_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v4y_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v5y_brick_dipole,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v0z_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v1z_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v2z_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v3z_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v4z_brick_spin,nzlo_out,nylo_out,nxlo_out); - memory->destroy3d_offset(v5z_brick_spin,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v0z_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v1z_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v2z_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v3z_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v4z_brick_dipole,nzlo_out,nylo_out,nxlo_out); + memory->destroy3d_offset(v5z_brick_dipole,nzlo_out,nylo_out,nxlo_out); - delete cg_peratom_spin; + delete cg_peratom_dipole; } /* ---------------------------------------------------------------------- - set global size of PPPMSpin grid = nx,ny,nz_pppm + set global size of PPPMDipole grid = nx,ny,nz_pppm used for charge accumulation, FFTs, and electric field interpolation ------------------------------------------------------------------------- */ -void PPPMSpin::set_grid_global() +void PPPMDipole::set_grid_global(double dipole2) { // use xprd,yprd,zprd - // adjust z dimension for 2d slab PPPMSpin - // 3d PPPMSpin just uses zprd since slab_volfactor = 1.0 + // adjust z dimension for 2d slab PPPMDipole + // 3d PPPMDipole just uses zprd since slab_volfactor = 1.0 double xprd = domain->xprd; double yprd = domain->yprd; @@ -817,14 +813,16 @@ void PPPMSpin::set_grid_global() if (!gewaldflag) { if (accuracy <= 0.0) error->all(FLERR,"KSpace accuracy must be > 0"); - if (sp2 == 0.0) - error->all(FLERR,"Must use kspace_modify gewald for systems with no spins"); + //if (mu2 == 0.0) + if (dipole2 == 0.0) + error->all(FLERR,"Must use kspace_modify gewald for systems with no dipoles"); g_ewald = (1.35 - 0.15*log(accuracy))/cutoff; //Try Newton Solver double g_ewald_new = - find_gewald_spin(g_ewald,cutoff,natoms,xprd*yprd*zprd,sp2); + find_gewald_dipole(g_ewald,cutoff,natoms,xprd*yprd*zprd,dipole2); + //find_gewald_dipole(g_ewald,cutoff,natoms,xprd*yprd*zprd,mu2); if (g_ewald_new > 0.0) g_ewald = g_ewald_new; - else error->warning(FLERR,"PPPMSpin spin Newton solver failed, " + else error->warning(FLERR,"PPPMDipole dipole Newton solver failed, " "using old method to estimate g_ewald"); } @@ -864,7 +862,7 @@ void PPPMSpin::set_grid_global() nzlo_fft = me_z*nz_pppm/npez_fft; nzhi_fft = (me_z+1)*nz_pppm/npez_fft - 1; - double df_kspace = compute_df_kspace_spin(); + double df_kspace = compute_df_kspace_dipole(dipole2); count++; @@ -889,31 +887,32 @@ void PPPMSpin::set_grid_global() h_z = zprd_slab/nz_pppm; if (nx_pppm >= OFFSET || ny_pppm >= OFFSET || nz_pppm >= OFFSET) - error->all(FLERR,"PPPMSpin grid is too large"); + error->all(FLERR,"PPPMDipole grid is too large"); } /* ---------------------------------------------------------------------- - compute estimated kspace force error for spins + compute estimated kspace force error for dipoles ------------------------------------------------------------------------- */ -double PPPMSpin::compute_df_kspace_spin() +double PPPMDipole::compute_df_kspace_dipole(double dipole2) { double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; double zprd_slab = zprd*slab_volfactor; bigint natoms = atom->natoms; - double qopt = compute_qopt_spin(); - double df_kspace = sqrt(qopt/natoms)*sp2/(3.0*xprd*yprd*zprd_slab); + double qopt = compute_qopt_dipole(); + //double df_kspace = sqrt(qopt/natoms)*mu2/(3.0*xprd*yprd*zprd_slab); + double df_kspace = sqrt(qopt/natoms)*dipole2/(3.0*xprd*yprd*zprd_slab); return df_kspace; } /* ---------------------------------------------------------------------- - compute qopt for spins with ik differentiation + compute qopt for dipoles with ik differentiation ------------------------------------------------------------------------- */ -double PPPMSpin::compute_qopt_spin() +double PPPMDipole::compute_qopt_dipole() { double qopt = 0.0; const double * const prd = domain->prd; @@ -984,7 +983,7 @@ double PPPMSpin::compute_qopt_spin() dot1 = unitkx*kper*qx + unitky*lper*qy + unitkz*mper*qz; dot2 = qx*qx + qy*qy + qz*qz; - //dot1 = dot1*dot1*dot1; // power of 3 for spin forces + //dot1 = dot1*dot1*dot1; // power of 3 for dipole forces //dot2 = dot2*dot2*dot2; u1 = sx*sy*sz; const double w2 = wx*wy*wz; @@ -1009,7 +1008,7 @@ double PPPMSpin::compute_qopt_spin() pre-compute modified (Hockney-Eastwood) Coulomb Green's function ------------------------------------------------------------------------- */ -void PPPMSpin::compute_gf_spin() +void PPPMDipole::compute_gf_dipole() { const double * const prd = domain->prd; @@ -1102,34 +1101,34 @@ void PPPMSpin::compute_gf_spin() calculate f(x) for use in Newton-Raphson solver ------------------------------------------------------------------------- */ -double PPPMSpin::newton_raphson_f() -{ - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - bigint natoms = atom->natoms; - - double df_rspace,df_kspace; - double vol = xprd*yprd*zprd; - double a = cutoff*g_ewald; - double rg2 = a*a; - double rg4 = rg2*rg2; - double rg6 = rg4*rg2; - double Cc = 4.0*rg4 + 6.0*rg2 + 3.0; - double Dc = 8.0*rg6 + 20.0*rg4 + 30.0*rg2 + 15.0; - df_rspace = (sp2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) * - sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) * exp(-rg2)); - df_kspace = compute_df_kspace_spin(); - - return df_rspace - df_kspace; -} +//double PPPMDipole::newton_raphson_f() +//{ +// double xprd = domain->xprd; +// double yprd = domain->yprd; +// double zprd = domain->zprd; +// bigint natoms = atom->natoms; +// +// double df_rspace,df_kspace; +// double vol = xprd*yprd*zprd; +// double a = cutoff*g_ewald; +// double rg2 = a*a; +// double rg4 = rg2*rg2; +// double rg6 = rg4*rg2; +// double Cc = 4.0*rg4 + 6.0*rg2 + 3.0; +// double Dc = 8.0*rg6 + 20.0*rg4 + 30.0*rg2 + 15.0; +// df_rspace = (mu2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) * +// sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) * exp(-rg2)); +// df_kspace = compute_df_kspace_dipole(); +// +// return df_rspace - df_kspace; +//} /* ---------------------------------------------------------------------- - find g_ewald parameter for spins based on desired accuracy + find g_ewald parameter for dipoles based on desired accuracy using a Newton-Raphson solver ------------------------------------------------------------------------- */ -double PPPMSpin::find_gewald_spin(double x, double Rc, +double PPPMDipole::find_gewald_dipole(double x, double Rc, bigint natoms, double vol, double b2) { double dx,tol; @@ -1141,7 +1140,7 @@ double PPPMSpin::find_gewald_spin(double x, double Rc, //Begin algorithm for (int i = 0; i < maxit; i++) { - dx = newton_raphson_f_spin(x,Rc,natoms,vol,b2) / derivf_spin(x,Rc,natoms,vol,b2); + dx = newton_raphson_f_dipole(x,Rc,natoms,vol,b2) / derivf_dipole(x,Rc,natoms,vol,b2); x = x - dx; //Update x if (fabs(dx) < tol) return x; if (x < 0 || x != x) // solver failed @@ -1151,10 +1150,10 @@ double PPPMSpin::find_gewald_spin(double x, double Rc, } /* ---------------------------------------------------------------------- - calculate f(x) objective function for spins + calculate f(x) objective function for dipoles ------------------------------------------------------------------------- */ -double PPPMSpin::newton_raphson_f_spin(double x, double Rc, bigint +double PPPMDipole::newton_raphson_f_dipole(double x, double Rc, bigint natoms, double vol, double b2) { double a = Rc*x; @@ -1171,21 +1170,21 @@ natoms, double vol, double b2) } /* ---------------------------------------------------------------------- - calculate numerical derivative f'(x) of objective function for spins + calculate numerical derivative f'(x) of objective function for dipoles ------------------------------------------------------------------------- */ -double PPPMSpin::derivf_spin(double x, double Rc, +double PPPMDipole::derivf_dipole(double x, double Rc, bigint natoms, double vol, double b2) { double h = 0.000001; //Derivative step-size - return (newton_raphson_f_spin(x + h,Rc,natoms,vol,b2) - newton_raphson_f_spin(x,Rc,natoms,vol,b2)) / h; + return (newton_raphson_f_dipole(x + h,Rc,natoms,vol,b2) - newton_raphson_f_dipole(x,Rc,natoms,vol,b2)) / h; } /* ---------------------------------------------------------------------- calculate the final estimate of the accuracy ------------------------------------------------------------------------- */ -double PPPMSpin::final_accuracy_spin() +double PPPMDipole::final_accuracy_dipole(double dipole2) { double xprd = domain->xprd; double yprd = domain->yprd; @@ -1194,7 +1193,7 @@ double PPPMSpin::final_accuracy_spin() bigint natoms = atom->natoms; if (natoms == 0) natoms = 1; // avoid division by zero - double df_kspace = compute_df_kspace_spin(); + double df_kspace = compute_df_kspace_dipole(mu2); double a = cutoff*g_ewald; double rg2 = a*a; @@ -1202,7 +1201,10 @@ double PPPMSpin::final_accuracy_spin() double rg6 = rg4*rg2; double Cc = 4.0*rg4 + 6.0*rg2 + 3.0; double Dc = 8.0*rg6 + 20.0*rg4 + 30.0*rg2 + 15.0; - double df_rspace = (sp2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) * + //double df_rspace = (mu2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) * + // sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) * + // exp(-rg2)); + double df_rspace = (dipole2/(sqrt(vol*powint(g_ewald,4)*powint(cutoff,9)*natoms)) * sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) * exp(-rg2)); @@ -1215,10 +1217,10 @@ double PPPMSpin::final_accuracy_spin() pre-compute Green's function denominator expansion coeffs, Gamma(2n) ------------------------------------------------------------------------- */ -void PPPMSpin::compute_gf_denom() +void PPPMDipole::compute_gf_denom() { if (gf_b) memory->destroy(gf_b); - memory->create(gf_b,order,"pppm_spin:gf_b"); + memory->create(gf_b,order,"pppm_dipole:gf_b"); int k,l,m; @@ -1244,7 +1246,7 @@ void PPPMSpin::compute_gf_denom() in global grid ------------------------------------------------------------------------- */ -void PPPMSpin::make_rho_spin() +void PPPMDipole::make_rho_dipole() { int l,m,n,nx,ny,nz,mx,my,mz; FFT_SCALAR dx,dy,dz; @@ -1254,11 +1256,11 @@ void PPPMSpin::make_rho_spin() // clear 3d density array - memset(&(densityx_brick_spin[nzlo_out][nylo_out][nxlo_out]),0, + memset(&(densityx_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0, ngrid*sizeof(FFT_SCALAR)); - memset(&(densityy_brick_spin[nzlo_out][nylo_out][nxlo_out]),0, + memset(&(densityy_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0, ngrid*sizeof(FFT_SCALAR)); - memset(&(densityz_brick_spin[nzlo_out][nylo_out][nxlo_out]),0, + memset(&(densityz_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0, ngrid*sizeof(FFT_SCALAR)); // loop over my charges, add their contribution to nearby grid points @@ -1266,8 +1268,7 @@ void PPPMSpin::make_rho_spin() // (dx,dy,dz) = distance to "lower left" grid pt // (mx,my,mz) = global coords of moving stencil pt - double **sp = atom->sp; - double spx,spy,spz; + double **mu = atom->mu; double **x = atom->x; int nlocal = atom->nlocal; @@ -1282,12 +1283,9 @@ void PPPMSpin::make_rho_spin() compute_rho1d(dx,dy,dz); - spx = sp[i][0]*sp[i][3]; - spy = sp[i][1]*sp[i][3]; - spz = sp[i][2]*sp[i][3]; - z0 = delvolinv * spx; - z1 = delvolinv * spy; - z2 = delvolinv * spz; + z0 = delvolinv * mu[i][0]; + z1 = delvolinv * mu[i][1]; + z2 = delvolinv * mu[i][2]; for (n = nlower; n <= nupper; n++) { mz = n+nz; y0 = z0*rho1d[2][n]; @@ -1300,9 +1298,9 @@ void PPPMSpin::make_rho_spin() x2 = y2*rho1d[1][m]; for (l = nlower; l <= nupper; l++) { mx = l+nx; - densityx_brick_spin[mz][my][mx] += x0*rho1d[0][l]; - densityy_brick_spin[mz][my][mx] += x1*rho1d[0][l]; - densityz_brick_spin[mz][my][mx] += x2*rho1d[0][l]; + densityx_brick_dipole[mz][my][mx] += x0*rho1d[0][l]; + densityy_brick_dipole[mz][my][mx] += x1*rho1d[0][l]; + densityz_brick_dipole[mz][my][mx] += x2*rho1d[0][l]; } } } @@ -1313,7 +1311,7 @@ void PPPMSpin::make_rho_spin() remap density from 3d brick decomposition to FFT decomposition ------------------------------------------------------------------------- */ -void PPPMSpin::brick2fft_spin() +void PPPMDipole::brick2fft_dipole() { int n,ix,iy,iz; @@ -1325,36 +1323,36 @@ void PPPMSpin::brick2fft_spin() for (iz = nzlo_in; iz <= nzhi_in; iz++) for (iy = nylo_in; iy <= nyhi_in; iy++) for (ix = nxlo_in; ix <= nxhi_in; ix++) { - densityx_fft_spin[n] = densityx_brick_spin[iz][iy][ix]; - densityy_fft_spin[n] = densityy_brick_spin[iz][iy][ix]; - densityz_fft_spin[n] = densityz_brick_spin[iz][iy][ix]; + densityx_fft_dipole[n] = densityx_brick_dipole[iz][iy][ix]; + densityy_fft_dipole[n] = densityy_brick_dipole[iz][iy][ix]; + densityz_fft_dipole[n] = densityz_brick_dipole[iz][iy][ix]; n++; } - remap->perform(densityx_fft_spin,densityx_fft_spin,work1); - remap->perform(densityy_fft_spin,densityy_fft_spin,work1); - remap->perform(densityz_fft_spin,densityz_fft_spin,work1); + remap->perform(densityx_fft_dipole,densityx_fft_dipole,work1); + remap->perform(densityy_fft_dipole,densityy_fft_dipole,work1); + remap->perform(densityz_fft_dipole,densityz_fft_dipole,work1); } /* ---------------------------------------------------------------------- FFT-based Poisson solver for ik ------------------------------------------------------------------------- */ -void PPPMSpin::poisson_ik_spin() +void PPPMDipole::poisson_ik_dipole() { int i,j,k,n,ii; double eng; double wreal,wimg; - // transform spin density (r -> k) + // transform dipole density (r -> k) n = 0; for (i = 0; i < nfft; i++) { - work1[n] = densityx_fft_spin[i]; + work1[n] = densityx_fft_dipole[i]; work1[n+1] = ZEROF; - work2[n] = densityy_fft_spin[i]; + work2[n] = densityy_fft_dipole[i]; work2[n+1] = ZEROF; - work3[n] = densityz_fft_spin[i]; + work3[n] = densityz_fft_dipole[i]; work3[n+1] = ZEROF; n += 2; } @@ -1421,7 +1419,7 @@ void PPPMSpin::poisson_ik_spin() // extra FFTs for per-atom energy/virial - if (vflag_atom) poisson_peratom_spin(); + if (vflag_atom) poisson_peratom_dipole(); // compute electric potential // FFT leaves data in 3d brick decomposition @@ -1443,7 +1441,7 @@ void PPPMSpin::poisson_ik_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - ux_brick_spin[k][j][i] = work4[n]; + ux_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1464,7 +1462,7 @@ void PPPMSpin::poisson_ik_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - uy_brick_spin[k][j][i] = work4[n]; + uy_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1485,7 +1483,7 @@ void PPPMSpin::poisson_ik_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - uz_brick_spin[k][j][i] = work4[n]; + uz_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1506,7 +1504,7 @@ void PPPMSpin::poisson_ik_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - vdxx_brick_spin[k][j][i] = work4[n]; + vdxx_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1527,7 +1525,7 @@ void PPPMSpin::poisson_ik_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - vdyy_brick_spin[k][j][i] = work4[n]; + vdyy_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1548,7 +1546,7 @@ void PPPMSpin::poisson_ik_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - vdzz_brick_spin[k][j][i] = work4[n]; + vdzz_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1569,7 +1567,7 @@ void PPPMSpin::poisson_ik_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - vdxy_brick_spin[k][j][i] = work4[n]; + vdxy_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1590,7 +1588,7 @@ void PPPMSpin::poisson_ik_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - vdxz_brick_spin[k][j][i] = work4[n]; + vdxz_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1611,7 +1609,7 @@ void PPPMSpin::poisson_ik_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - vdyz_brick_spin[k][j][i] = work4[n]; + vdyz_brick_dipole[k][j][i] = work4[n]; n += 2; } } @@ -1620,7 +1618,7 @@ void PPPMSpin::poisson_ik_spin() FFT-based Poisson solver for per-atom energy/virial ------------------------------------------------------------------------- */ -void PPPMSpin::poisson_peratom_spin() +void PPPMDipole::poisson_peratom_dipole() { int i,ii,j,k,n; @@ -1647,7 +1645,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v0x_brick_spin[k][j][i] = work4[n]; + v0x_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1670,7 +1668,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v0y_brick_spin[k][j][i] = work4[n]; + v0y_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1693,7 +1691,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v0z_brick_spin[k][j][i] = work4[n]; + v0z_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1716,7 +1714,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v1x_brick_spin[k][j][i] = work4[n]; + v1x_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1739,7 +1737,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v1y_brick_spin[k][j][i] = work4[n]; + v1y_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1762,7 +1760,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v1z_brick_spin[k][j][i] = work4[n]; + v1z_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1785,7 +1783,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v2x_brick_spin[k][j][i] = work4[n]; + v2x_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1808,7 +1806,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v2y_brick_spin[k][j][i] = work4[n]; + v2y_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1831,7 +1829,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v2z_brick_spin[k][j][i] = work4[n]; + v2z_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1854,7 +1852,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v3x_brick_spin[k][j][i] = work4[n]; + v3x_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1877,7 +1875,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v3y_brick_spin[k][j][i] = work4[n]; + v3y_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1900,7 +1898,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v3z_brick_spin[k][j][i] = work4[n]; + v3z_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1923,7 +1921,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v4x_brick_spin[k][j][i] = work4[n]; + v4x_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1946,7 +1944,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v4y_brick_spin[k][j][i] = work4[n]; + v4y_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1969,7 +1967,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v4z_brick_spin[k][j][i] = work4[n]; + v4z_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -1992,7 +1990,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v5x_brick_spin[k][j][i] = work4[n]; + v5x_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -2015,7 +2013,7 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v5y_brick_spin[k][j][i] = work4[n]; + v5y_brick_dipole[k][j][i] = work4[n]; n += 2; } @@ -2038,16 +2036,16 @@ void PPPMSpin::poisson_peratom_spin() for (k = nzlo_in; k <= nzhi_in; k++) for (j = nylo_in; j <= nyhi_in; j++) for (i = nxlo_in; i <= nxhi_in; i++) { - v5z_brick_spin[k][j][i] = work4[n]; + v5z_brick_dipole[k][j][i] = work4[n]; n += 2; } } /* ---------------------------------------------------------------------- - interpolate from grid to get magnetic field & force on my particles for ik + interpolate from grid to get electric field & force on my particles for ik ------------------------------------------------------------------------- */ -void PPPMSpin::fieldforce_ik_spin() +void PPPMDipole::fieldforce_ik_dipole() { int i,l,m,n,nx,ny,nz,mx,my,mz; FFT_SCALAR dx,dy,dz; @@ -2060,11 +2058,11 @@ void PPPMSpin::fieldforce_ik_spin() // (dx,dy,dz) = distance to "lower left" grid pt // (mx,my,mz) = global coords of moving stencil pt - double **sp = atom->sp; - double spx,spy,spz; + + double **mu = atom->mu; double **x = atom->x; double **f = atom->f; - double **fm = atom->fm; + double **t = atom->torque; int nlocal = atom->nlocal; @@ -2089,36 +2087,29 @@ void PPPMSpin::fieldforce_ik_spin() for (l = nlower; l <= nupper; l++) { mx = l+nx; x0 = y0*rho1d[0][l]; - ex -= x0*ux_brick_spin[mz][my][mx]; - ey -= x0*uy_brick_spin[mz][my][mx]; - ez -= x0*uz_brick_spin[mz][my][mx]; - vxx -= x0*vdxx_brick_spin[mz][my][mx]; - vyy -= x0*vdyy_brick_spin[mz][my][mx]; - vzz -= x0*vdzz_brick_spin[mz][my][mx]; - vxy -= x0*vdxy_brick_spin[mz][my][mx]; - vxz -= x0*vdxz_brick_spin[mz][my][mx]; - vyz -= x0*vdyz_brick_spin[mz][my][mx]; + ex -= x0*ux_brick_dipole[mz][my][mx]; + ey -= x0*uy_brick_dipole[mz][my][mx]; + ez -= x0*uz_brick_dipole[mz][my][mx]; + vxx -= x0*vdxx_brick_dipole[mz][my][mx]; + vyy -= x0*vdyy_brick_dipole[mz][my][mx]; + vzz -= x0*vdzz_brick_dipole[mz][my][mx]; + vxy -= x0*vdxy_brick_dipole[mz][my][mx]; + vxz -= x0*vdxz_brick_dipole[mz][my][mx]; + vyz -= x0*vdyz_brick_dipole[mz][my][mx]; } } } - // convert M-field to mech. and mag. forces + // convert E-field to torque - const double spfactor = mub2mu0 * scale; - spx = sp[i][0]*sp[i][3]; - spy = sp[i][1]*sp[i][3]; - spz = sp[i][2]*sp[i][3]; - f[i][0] += spfactor*(vxx*spx + vxy*spy + vxz*spz); - f[i][1] += spfactor*(vxy*spx + vyy*spy + vyz*spz); - f[i][2] += spfactor*(vxz*spx + vyz*spy + vzz*spz); - - const double spfactorh = mub2mu0hbinv * scale; - fm[i][0] += spfactorh*ex; - fm[i][1] += spfactorh*ey; - fm[i][2] += spfactorh*ez; - - // create a new vector (in atom_spin style ?) to store long-range fm tables + const double mufactor = qqrd2e * scale; + f[i][0] += mufactor*(vxx*mu[i][0] + vxy*mu[i][1] + vxz*mu[i][2]); + f[i][1] += mufactor*(vxy*mu[i][0] + vyy*mu[i][1] + vyz*mu[i][2]); + f[i][2] += mufactor*(vxz*mu[i][0] + vyz*mu[i][1] + vzz*mu[i][2]); + t[i][0] += mufactor*(mu[i][1]*ez - mu[i][2]*ey); + t[i][1] += mufactor*(mu[i][2]*ex - mu[i][0]*ez); + t[i][2] += mufactor*(mu[i][0]*ey - mu[i][1]*ex); } } @@ -2126,7 +2117,7 @@ void PPPMSpin::fieldforce_ik_spin() interpolate from grid to get per-atom energy/virial ------------------------------------------------------------------------- */ -void PPPMSpin::fieldforce_peratom_spin() +void PPPMDipole::fieldforce_peratom_dipole() { int i,l,m,n,nx,ny,nz,mx,my,mz; FFT_SCALAR dx,dy,dz,x0,y0,z0; @@ -2140,8 +2131,7 @@ void PPPMSpin::fieldforce_peratom_spin() // (dx,dy,dz) = distance to "lower left" grid pt // (mx,my,mz) = global coords of moving stencil pt - double **sp = atom->sp; - double spx,spy,spz; + double **mu = atom->mu; double **x = atom->x; int nlocal = atom->nlocal; @@ -2170,45 +2160,42 @@ void PPPMSpin::fieldforce_peratom_spin() mx = l+nx; x0 = y0*rho1d[0][l]; if (eflag_atom) { - ux += x0*ux_brick_spin[mz][my][mx]; - uy += x0*uy_brick_spin[mz][my][mx]; - uz += x0*uz_brick_spin[mz][my][mx]; + ux += x0*ux_brick_dipole[mz][my][mx]; + uy += x0*uy_brick_dipole[mz][my][mx]; + uz += x0*uz_brick_dipole[mz][my][mx]; } if (vflag_atom) { - v0x += x0*v0x_brick_spin[mz][my][mx]; - v1x += x0*v1x_brick_spin[mz][my][mx]; - v2x += x0*v2x_brick_spin[mz][my][mx]; - v3x += x0*v3x_brick_spin[mz][my][mx]; - v4x += x0*v4x_brick_spin[mz][my][mx]; - v5x += x0*v5x_brick_spin[mz][my][mx]; - v0y += x0*v0y_brick_spin[mz][my][mx]; - v1y += x0*v1y_brick_spin[mz][my][mx]; - v2y += x0*v2y_brick_spin[mz][my][mx]; - v3y += x0*v3y_brick_spin[mz][my][mx]; - v4y += x0*v4y_brick_spin[mz][my][mx]; - v5y += x0*v5y_brick_spin[mz][my][mx]; - v0z += x0*v0z_brick_spin[mz][my][mx]; - v1z += x0*v1z_brick_spin[mz][my][mx]; - v2z += x0*v2z_brick_spin[mz][my][mx]; - v3z += x0*v3z_brick_spin[mz][my][mx]; - v4z += x0*v4z_brick_spin[mz][my][mx]; - v5z += x0*v5z_brick_spin[mz][my][mx]; + v0x += x0*v0x_brick_dipole[mz][my][mx]; + v1x += x0*v1x_brick_dipole[mz][my][mx]; + v2x += x0*v2x_brick_dipole[mz][my][mx]; + v3x += x0*v3x_brick_dipole[mz][my][mx]; + v4x += x0*v4x_brick_dipole[mz][my][mx]; + v5x += x0*v5x_brick_dipole[mz][my][mx]; + v0y += x0*v0y_brick_dipole[mz][my][mx]; + v1y += x0*v1y_brick_dipole[mz][my][mx]; + v2y += x0*v2y_brick_dipole[mz][my][mx]; + v3y += x0*v3y_brick_dipole[mz][my][mx]; + v4y += x0*v4y_brick_dipole[mz][my][mx]; + v5y += x0*v5y_brick_dipole[mz][my][mx]; + v0z += x0*v0z_brick_dipole[mz][my][mx]; + v1z += x0*v1z_brick_dipole[mz][my][mx]; + v2z += x0*v2z_brick_dipole[mz][my][mx]; + v3z += x0*v3z_brick_dipole[mz][my][mx]; + v4z += x0*v4z_brick_dipole[mz][my][mx]; + v5z += x0*v5z_brick_dipole[mz][my][mx]; } } } } - spx = sp[i][0]*sp[i][3]; - spy = sp[i][1]*sp[i][3]; - spz = sp[i][2]*sp[i][3]; - if (eflag_atom) eatom[i] += spx*ux + spy*uy + spz*uz; + if (eflag_atom) eatom[i] += mu[i][0]*ux + mu[i][1]*uy + mu[i][2]*uz; if (vflag_atom) { - vatom[i][0] += spx*v0x + spy*v0y + spz*v0z; - vatom[i][1] += spx*v1x + spy*v1y + spz*v1z; - vatom[i][2] += spx*v2x + spy*v2y + spz*v2z; - vatom[i][3] += spx*v3x + spy*v3y + spz*v3z; - vatom[i][4] += spx*v4x + spy*v4y + spz*v4z; - vatom[i][5] += spx*v5x + spy*v5y + spz*v5z; + vatom[i][0] += mu[i][0]*v0x + mu[i][1]*v0y + mu[i][2]*v0z; + vatom[i][1] += mu[i][0]*v1x + mu[i][1]*v1y + mu[i][2]*v1z; + vatom[i][2] += mu[i][0]*v2x + mu[i][1]*v2y + mu[i][2]*v2z; + vatom[i][3] += mu[i][0]*v3x + mu[i][1]*v3y + mu[i][2]*v3z; + vatom[i][4] += mu[i][0]*v4x + mu[i][1]*v4y + mu[i][2]*v4z; + vatom[i][5] += mu[i][0]*v5x + mu[i][1]*v5y + mu[i][2]*v5z; } } } @@ -2217,20 +2204,20 @@ void PPPMSpin::fieldforce_peratom_spin() pack own values to buf to send to another proc ------------------------------------------------------------------------- */ -void PPPMSpin::pack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list) +void PPPMDipole::pack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list) { int n = 0; - if (flag == FORWARD_SP) { - FFT_SCALAR *src_ux = &ux_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *src_uy = &uy_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *src_uz = &uz_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *src_vxx = &vdxx_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *src_vyy = &vdyy_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *src_vzz = &vdzz_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *src_vxy = &vdxy_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *src_vxz = &vdxz_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *src_vyz = &vdyz_brick_spin[nzlo_out][nylo_out][nxlo_out]; + if (flag == FORWARD_MU) { + FFT_SCALAR *src_ux = &ux_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *src_uy = &uy_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *src_uz = &uz_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *src_vxx = &vdxx_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *src_vyy = &vdyy_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *src_vzz = &vdzz_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *src_vxy = &vdxy_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *src_vxz = &vdxz_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *src_vyz = &vdyz_brick_dipole[nzlo_out][nylo_out][nxlo_out]; for (int i = 0; i < nlist; i++) { buf[n++] = src_ux[list[i]]; buf[n++] = src_uy[list[i]]; @@ -2242,25 +2229,25 @@ void PPPMSpin::pack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list) buf[n++] = src_vxz[list[i]]; buf[n++] = src_vyz[list[i]]; } - } else if (flag == FORWARD_SP_PERATOM) { - FFT_SCALAR *v0xsrc = &v0x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v1xsrc = &v1x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v2xsrc = &v2x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v3xsrc = &v3x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v4xsrc = &v4x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v5xsrc = &v5x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v0ysrc = &v0y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v1ysrc = &v1y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v2ysrc = &v2y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v3ysrc = &v3y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v4ysrc = &v4y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v5ysrc = &v5y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v0zsrc = &v0z_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v1zsrc = &v1z_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v2zsrc = &v2z_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v3zsrc = &v3z_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v4zsrc = &v4z_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v5zsrc = &v5z_brick_spin[nzlo_out][nylo_out][nxlo_out]; + } else if (flag == FORWARD_MU_PERATOM) { + FFT_SCALAR *v0xsrc = &v0x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v1xsrc = &v1x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v2xsrc = &v2x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v3xsrc = &v3x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v4xsrc = &v4x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v5xsrc = &v5x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v0ysrc = &v0y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v1ysrc = &v1y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v2ysrc = &v2y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v3ysrc = &v3y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v4ysrc = &v4y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v5ysrc = &v5y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v0zsrc = &v0z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v1zsrc = &v1z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v2zsrc = &v2z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v3zsrc = &v3z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v4zsrc = &v4z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v5zsrc = &v5z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; for (int i = 0; i < nlist; i++) { buf[n++] = v0xsrc[list[i]]; buf[n++] = v1xsrc[list[i]]; @@ -2288,20 +2275,20 @@ void PPPMSpin::pack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list) unpack another proc's own values from buf and set own ghost values ------------------------------------------------------------------------- */ -void PPPMSpin::unpack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list) +void PPPMDipole::unpack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list) { int n = 0; - if (flag == FORWARD_SP) { - FFT_SCALAR *dest_ux = &ux_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *dest_uy = &uy_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *dest_uz = &uz_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *dest_vxx = &vdxx_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *dest_vyy = &vdyy_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *dest_vzz = &vdzz_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *dest_vxy = &vdxy_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *dest_vxz = &vdxz_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *dest_vyz = &vdyz_brick_spin[nzlo_out][nylo_out][nxlo_out]; + if (flag == FORWARD_MU) { + FFT_SCALAR *dest_ux = &ux_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *dest_uy = &uy_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *dest_uz = &uz_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *dest_vxx = &vdxx_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *dest_vyy = &vdyy_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *dest_vzz = &vdzz_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *dest_vxy = &vdxy_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *dest_vxz = &vdxz_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *dest_vyz = &vdyz_brick_dipole[nzlo_out][nylo_out][nxlo_out]; for (int i = 0; i < nlist; i++) { dest_ux[list[i]] = buf[n++]; dest_uy[list[i]] = buf[n++]; @@ -2313,25 +2300,25 @@ void PPPMSpin::unpack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list) dest_vxz[list[i]] = buf[n++]; dest_vyz[list[i]] = buf[n++]; } - } else if (flag == FORWARD_SP_PERATOM) { - FFT_SCALAR *v0xsrc = &v0x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v1xsrc = &v1x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v2xsrc = &v2x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v3xsrc = &v3x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v4xsrc = &v4x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v5xsrc = &v5x_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v0ysrc = &v0y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v1ysrc = &v1y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v2ysrc = &v2y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v3ysrc = &v3y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v4ysrc = &v4y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v5ysrc = &v5y_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v0zsrc = &v0z_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v1zsrc = &v1z_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v2zsrc = &v2z_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v3zsrc = &v3z_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v4zsrc = &v4z_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *v5zsrc = &v5z_brick_spin[nzlo_out][nylo_out][nxlo_out]; + } else if (flag == FORWARD_MU_PERATOM) { + FFT_SCALAR *v0xsrc = &v0x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v1xsrc = &v1x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v2xsrc = &v2x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v3xsrc = &v3x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v4xsrc = &v4x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v5xsrc = &v5x_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v0ysrc = &v0y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v1ysrc = &v1y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v2ysrc = &v2y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v3ysrc = &v3y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v4ysrc = &v4y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v5ysrc = &v5y_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v0zsrc = &v0z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v1zsrc = &v1z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v2zsrc = &v2z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v3zsrc = &v3z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v4zsrc = &v4z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *v5zsrc = &v5z_brick_dipole[nzlo_out][nylo_out][nxlo_out]; for (int i = 0; i < nlist; i++) { v0xsrc[list[i]] = buf[n++]; v1xsrc[list[i]] = buf[n++]; @@ -2359,17 +2346,17 @@ void PPPMSpin::unpack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list) pack ghost values into buf to send to another proc ------------------------------------------------------------------------- */ -void PPPMSpin::pack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list) +void PPPMDipole::pack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list) { int n = 0; - if (flag == REVERSE_SP) { - FFT_SCALAR *src_spin0 = &densityx_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *src_spin1 = &densityy_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *src_spin2 = &densityz_brick_spin[nzlo_out][nylo_out][nxlo_out]; + if (flag == REVERSE_MU) { + FFT_SCALAR *src_dipole0 = &densityx_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *src_dipole1 = &densityy_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *src_dipole2 = &densityz_brick_dipole[nzlo_out][nylo_out][nxlo_out]; for (int i = 0; i < nlist; i++) { - buf[n++] = src_spin0[list[i]]; - buf[n++] = src_spin1[list[i]]; - buf[n++] = src_spin2[list[i]]; + buf[n++] = src_dipole0[list[i]]; + buf[n++] = src_dipole1[list[i]]; + buf[n++] = src_dipole2[list[i]]; } } } @@ -2378,17 +2365,17 @@ void PPPMSpin::pack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list) unpack another proc's ghost values from buf and add to own values ------------------------------------------------------------------------- */ -void PPPMSpin::unpack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list) +void PPPMDipole::unpack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list) { int n = 0; - if (flag == REVERSE_SP) { - FFT_SCALAR *dest_spin0 = &densityx_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *dest_spin1 = &densityy_brick_spin[nzlo_out][nylo_out][nxlo_out]; - FFT_SCALAR *dest_spin2 = &densityz_brick_spin[nzlo_out][nylo_out][nxlo_out]; + if (flag == REVERSE_MU) { + FFT_SCALAR *dest_dipole0 = &densityx_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *dest_dipole1 = &densityy_brick_dipole[nzlo_out][nylo_out][nxlo_out]; + FFT_SCALAR *dest_dipole2 = &densityz_brick_dipole[nzlo_out][nylo_out][nxlo_out]; for (int i = 0; i < nlist; i++) { - dest_spin0[list[i]] += buf[n++]; - dest_spin1[list[i]] += buf[n++]; - dest_spin2[list[i]] += buf[n++]; + dest_dipole0[list[i]] += buf[n++]; + dest_dipole1[list[i]] += buf[n++]; + dest_dipole2[list[i]] += buf[n++]; } } } @@ -2401,50 +2388,57 @@ void PPPMSpin::unpack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list) extended to non-neutral systems (J. Chem. Phys. 131, 094107). ------------------------------------------------------------------------- */ -void PPPMSpin::slabcorr() +void PPPMDipole::slabcorr() { - // compute local contribution to global spin moment + // compute local contribution to global dipole moment double **x = atom->x; double zprd = domain->zprd; int nlocal = atom->nlocal; - double spin = 0.0; - double **sp = atom->sp; - double spx,spy,spz; - for (int i = 0; i < nlocal; i++) { - spz = sp[i][2]*sp[i][3]; - spin += spz; + double dipole = 0.0; + double **mu = atom->mu; + for (int i = 0; i < nlocal; i++) dipole += mu[i][2]; + + // sum local contributions to get global dipole moment + + double dipole_all; + MPI_Allreduce(&dipole,&dipole_all,1,MPI_DOUBLE,MPI_SUM,world); + + // need to make non-neutral systems and/or + // per-atom energy translationally invariant + + if (eflag_atom || fabs(qsum) > SMALL) { + + error->all(FLERR,"Cannot (yet) use kspace slab correction with " + "long-range dipoles and non-neutral systems or per-atom energy"); } - // sum local contributions to get global spin moment - - double spin_all; - MPI_Allreduce(&spin,&spin_all,1,MPI_DOUBLE,MPI_SUM,world); - // compute corrections - const double e_slabcorr = MY_2PI*(spin_all*spin_all/12.0)/volume; - const double spscale = mub2mu0 * scale; + const double e_slabcorr = MY_2PI*(dipole_all*dipole_all/12.0)/volume; + const double qscale = qqrd2e * scale; - if (eflag_global) energy += spscale * e_slabcorr; + if (eflag_global) energy += qscale * e_slabcorr; // per-atom energy if (eflag_atom) { - double efact = spscale * MY_2PI/volume/12.0; - for (int i = 0; i < nlocal; i++) { - spz = sp[i][2]*sp[i][3]; - eatom[i] += efact * spz * spin_all; - } + double efact = qscale * MY_2PI/volume/12.0; + for (int i = 0; i < nlocal; i++) + eatom[i] += efact * mu[i][2]*dipole_all; } - // add on mag. force corrections + // add on torque corrections - double ffact = spscale * (-4.0*MY_PI/volume); - double **fm = atom->fm; - for (int i = 0; i < nlocal; i++) { - fm[i][2] += ffact * spin_all; + if (atom->torque) { + double ffact = qscale * (-4.0*MY_PI/volume); + double **mu = atom->mu; + double **torque = atom->torque; + for (int i = 0; i < nlocal; i++) { + torque[i][0] += ffact * dipole_all * mu[i][1]; + torque[i][1] += -ffact * dipole_all * mu[i][0]; + } } } @@ -2452,7 +2446,7 @@ void PPPMSpin::slabcorr() perform and time the 1d FFTs required for N timesteps ------------------------------------------------------------------------- */ -int PPPMSpin::timing_1d(int n, double &time1d) +int PPPMDipole::timing_1d(int n, double &time1d) { double time1,time2; @@ -2487,7 +2481,7 @@ int PPPMSpin::timing_1d(int n, double &time1d) perform and time the 3d FFTs required for N timesteps ------------------------------------------------------------------------- */ -int PPPMSpin::timing_3d(int n, double &time3d) +int PPPMDipole::timing_3d(int n, double &time3d) { double time1,time2; @@ -2522,7 +2516,7 @@ int PPPMSpin::timing_3d(int n, double &time3d) memory usage of local arrays ------------------------------------------------------------------------- */ -double PPPMSpin::memory_usage() +double PPPMDipole::memory_usage() { double bytes = nmax*3 * sizeof(double); int nbrick = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) * @@ -2536,43 +2530,37 @@ double PPPMSpin::memory_usage() if (peratom_allocate_flag) bytes += 21 * nbrick * sizeof(FFT_SCALAR); - if (cg_spin) bytes += cg_spin->memory_usage(); - if (cg_peratom_spin) bytes += cg_peratom_spin->memory_usage(); + if (cg_dipole) bytes += cg_dipole->memory_usage(); + if (cg_peratom_dipole) bytes += cg_peratom_dipole->memory_usage(); return bytes; } /* ---------------------------------------------------------------------- - compute spsum,spsqsum,sp2 - called initially, when particle count changes, when spins are changed + compute musum,musqsum,mu2 + called initially, when particle count changes, when dipoles are changed ------------------------------------------------------------------------- */ -void PPPMSpin::spsum_spsq() +void PPPMDipole::musum_musq() { const int nlocal = atom->nlocal; - spsum = spsqsum = sp2 = 0.0; - if (atom->sp_flag) { - double **sp = atom->sp; - double spx, spy, spz; - double spsum_local(0.0), spsqsum_local(0.0); - - // not exactly the good loop: need to add norm of spins + musum = musqsum = mu2 = 0.0; + if (atom->mu_flag) { + double** mu = atom->mu; + double musum_local(0.0), musqsum_local(0.0); for (int i = 0; i < nlocal; i++) { - spx = sp[i][0]*sp[i][3]; - spy = sp[i][1]*sp[i][3]; - spz = sp[i][2]*sp[i][3]; - spsum_local += spx + spy + spz; - spsqsum_local += spx*spx + spy*spy + spz*spz; + musum_local += mu[i][0] + mu[i][1] + mu[i][2]; + musqsum_local += mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]; } - MPI_Allreduce(&spsum_local,&spsum,1,MPI_DOUBLE,MPI_SUM,world); - MPI_Allreduce(&spsqsum_local,&spsqsum,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&musum_local,&musum,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&musqsum_local,&musqsum,1,MPI_DOUBLE,MPI_SUM,world); - sp2 = spsqsum * mub2mu0; + mu2 = musqsum * force->qqrd2e; } - if (sp2 == 0 && comm->me == 0) - error->all(FLERR,"Using kspace solver PPPMSpin on system with no spins"); + if (mu2 == 0 && comm->me == 0) + error->all(FLERR,"Using kspace solver PPPMDipole on system with no dipoles"); } diff --git a/src/KSPACE/pppm_dipole.h b/src/KSPACE/pppm_dipole.h new file mode 100644 index 0000000000..8db28b540a --- /dev/null +++ b/src/KSPACE/pppm_dipole.h @@ -0,0 +1,213 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef KSPACE_CLASS + +KSpaceStyle(pppm/dipole,PPPMDipole) + +#else + +#ifndef LMP_PPPM_DIPOLE_H +#define LMP_PPPM_DIPOLE_H + +#include "pppm.h" + +namespace LAMMPS_NS { + +class PPPMDipole : public PPPM { + public: + PPPMDipole(class LAMMPS *, int, char **); + virtual ~PPPMDipole(); + void init(); + void setup(); + void setup_grid(); + void compute(int, int); + int timing_1d(int, double &); + int timing_3d(int, double &); + double memory_usage(); + + protected: + void set_grid_global(double); + //double newton_raphson_f(); + + void allocate(); + void allocate_peratom(); + void deallocate(); + void deallocate_peratom(); + void compute_gf_denom(); + + void slabcorr(); + + // grid communication + + void pack_forward(int, FFT_SCALAR *, int, int *); + void unpack_forward(int, FFT_SCALAR *, int, int *); + void pack_reverse(int, FFT_SCALAR *, int, int *); + void unpack_reverse(int, FFT_SCALAR *, int, int *); + + // dipole + + FFT_SCALAR ***densityx_brick_dipole,***densityy_brick_dipole,***densityz_brick_dipole; + FFT_SCALAR ***vdxx_brick_dipole,***vdyy_brick_dipole,***vdzz_brick_dipole; + FFT_SCALAR ***vdxy_brick_dipole,***vdxz_brick_dipole,***vdyz_brick_dipole; + FFT_SCALAR ***ux_brick_dipole,***uy_brick_dipole,***uz_brick_dipole; + FFT_SCALAR ***v0x_brick_dipole,***v1x_brick_dipole,***v2x_brick_dipole; + FFT_SCALAR ***v3x_brick_dipole,***v4x_brick_dipole,***v5x_brick_dipole; + FFT_SCALAR ***v0y_brick_dipole,***v1y_brick_dipole,***v2y_brick_dipole; + FFT_SCALAR ***v3y_brick_dipole,***v4y_brick_dipole,***v5y_brick_dipole; + FFT_SCALAR ***v0z_brick_dipole,***v1z_brick_dipole,***v2z_brick_dipole; + FFT_SCALAR ***v3z_brick_dipole,***v4z_brick_dipole,***v5z_brick_dipole; + FFT_SCALAR *work3,*work4; + FFT_SCALAR *densityx_fft_dipole,*densityy_fft_dipole,*densityz_fft_dipole; + class GridComm *cg_dipole; + class GridComm *cg_peratom_dipole; + int only_dipole_flag; + double musum,musqsum,mu2; + double find_gewald_dipole(double, double, bigint, double, double); + double newton_raphson_f_dipole(double, double, bigint, double, double); + double derivf_dipole(double, double, bigint, double, double); + double compute_df_kspace_dipole(double); + double compute_qopt_dipole(); + void compute_gf_dipole(); + void make_rho_dipole(); + void brick2fft_dipole(); + void poisson_ik_dipole(); + void poisson_peratom_dipole(); + void fieldforce_ik_dipole(); + void fieldforce_peratom_dipole(); + double final_accuracy_dipole(double dipole2); + void musum_musq(); + +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Cannot (yet) use charges with Kspace style PPPMDipole + +Charge-dipole interactions are not yet implemented in PPPMDipole so this +feature is not yet supported. + +E: Must redefine kspace_style after changing to triclinic box + +Self-explanatory. + +E: Kspace style requires atom attribute mu + +The atom style defined does not have this attribute. + +E: Cannot (yet) use kspace_modify diff ad with dipoles + +This feature is not yet supported. + +E: Cannot (yet) use 'electron' units with dipoles + +This feature is not yet supported. + +E: Cannot yet use triclinic cells with PPPMDipole + +This feature is not yet supported. + +E: Cannot yet use TIP4P with PPPMDipole + +This feature is not yet supported. + +E: Cannot use nonperiodic boundaries with PPPM + +For kspace style pppm, all 3 dimensions must have periodic boundaries +unless you use the kspace_modify command to define a 2d slab with a +non-periodic z dimension. + +E: Incorrect boundaries with slab PPPM + +Must have periodic x,y dimensions and non-periodic z dimension to use +2d slab option with PPPM. + +E: PPPM order cannot be < 2 or > than %d + +This is a limitation of the PPPM implementation in LAMMPS. + +E: KSpace style is incompatible with Pair style + +Setting a kspace style requires that a pair style with matching +long-range dipole components be used. + +W: Reducing PPPM order b/c stencil extends beyond nearest neighbor processor + +This may lead to a larger grid than desired. See the kspace_modify overlap +command to prevent changing of the PPPM order. + +E: PPPM order < minimum allowed order + +The default minimum order is 2. This can be reset by the +kspace_modify minorder command. + +E: PPPM grid stencil extends beyond nearest neighbor processor + +This is not allowed if the kspace_modify overlap setting is no. + +E: KSpace accuracy must be > 0 + +The kspace accuracy designated in the input must be greater than zero. + +E: Could not compute grid size + +The code is unable to compute a grid size consistent with the desired +accuracy. This error should not occur for typical problems. Please +send an email to the developers. + +E: PPPM grid is too large + +The global PPPM grid is larger than OFFSET in one or more dimensions. +OFFSET is currently set to 4096. You likely need to decrease the +requested accuracy. + +E: Could not compute g_ewald + +The Newton-Raphson solver failed to converge to a good value for +g_ewald. This error should not occur for typical problems. Please +send an email to the developers. + +E: Non-numeric box dimensions - simulation unstable + +The box size has apparently blown up. + +E: Out of range atoms - cannot compute PPPM + +One or more atoms are attempting to map their charge to a PPPM grid +point that is not owned by a processor. This is likely for one of two +reasons, both of them bad. First, it may mean that an atom near the +boundary of a processor's sub-domain has moved more than 1/2 the +"neighbor skin distance"_neighbor.html without neighbor lists being +rebuilt and atoms being migrated to new processors. This also means +you may be missing pairwise interactions that need to be computed. +The solution is to change the re-neighboring criteria via the +"neigh_modify"_neigh_modify command. The safest settings are "delay 0 +every 1 check yes". Second, it may mean that an atom has moved far +outside a processor's sub-domain or even the entire simulation box. +This indicates bad physics, e.g. due to highly overlapping atoms, too +large a timestep, etc. + +E: Using kspace solver PPPMDipole on system with no dipoles + +Must have non-zero dipoles with PPPMDipole. + +E: Must use kspace_modify gewald for system with no dipoles + +Self-explanatory. + +*/ diff --git a/src/KSPACE/pppm_dipole_spin.cpp b/src/KSPACE/pppm_dipole_spin.cpp new file mode 100644 index 0000000000..4fde7ba101 --- /dev/null +++ b/src/KSPACE/pppm_dipole_spin.cpp @@ -0,0 +1,750 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Stan Moore (SNL) + Julien Tranchida (SNL) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include +#include "pppm_dipole_spin.h" +#include "atom.h" +#include "comm.h" +#include "gridcomm.h" +#include "neighbor.h" +#include "force.h" +#include "pair.h" +#include "bond.h" +#include "angle.h" +#include "domain.h" +#include "fft3d_wrap.h" +#include "remap_wrap.h" +#include "memory.h" +#include "error.h" +#include "update.h" + +#include "math_const.h" +#include "math_special.h" + +using namespace LAMMPS_NS; +using namespace MathConst; +using namespace MathSpecial; + +#define MAXORDER 7 +#define OFFSET 16384 +#define LARGE 10000.0 +#define SMALL 0.00001 +#define EPS_HOC 1.0e-7 + +enum{REVERSE_SP}; +enum{FORWARD_SP,FORWARD_SP_PERATOM}; + +#ifdef FFT_SINGLE +#define ZEROF 0.0f +#define ONEF 1.0f +#else +#define ZEROF 0.0 +#define ONEF 1.0 +#endif + +/* ---------------------------------------------------------------------- */ + +PPPMDipoleSpin::PPPMDipoleSpin(LAMMPS *lmp, int narg, char **arg) : + PPPMDipole(lmp, narg, arg) +{ + dipoleflag = 0; + spinflag = 1; + group_group_enable = 0; + + cg_dipole = NULL; + cg_peratom_dipole = NULL; +} + +/* ---------------------------------------------------------------------- + free all memory +------------------------------------------------------------------------- */ + +PPPMDipoleSpin::~PPPMDipoleSpin() +{ + if (copymode) return; + + deallocate(); + if (peratom_allocate_flag) deallocate_peratom(); + fft1 = NULL; + fft2 = NULL; + remap = NULL; + cg_dipole = NULL; +} + +/* ---------------------------------------------------------------------- + called once before run +------------------------------------------------------------------------- */ + +void PPPMDipoleSpin::init() +{ + if (me == 0) { + if (screen) fprintf(screen,"PPPMDipoleSpin initialization ...\n"); + if (logfile) fprintf(logfile,"PPPMDipoleSpin initialization ...\n"); + } + + // error check + + spinflag = atom->sp?1:0; + + triclinic_check(); + + if (triclinic != domain->triclinic) + error->all(FLERR,"Must redefine kspace_style after changing to triclinic box"); + + if (domain->dimension == 2) error->all(FLERR, + "Cannot use PPPMDipoleSpin with 2d simulation"); + if (comm->style != 0) + error->universe_all(FLERR,"PPPMDipoleSpin can only currently be used with " + "comm_style brick"); + + if (!atom->sp) error->all(FLERR,"Kspace style requires atom attribute sp"); + + if (atom->sp && differentiation_flag == 1) error->all(FLERR,"Cannot (yet) use kspace_modify diff" + " ad with spins"); + + if (spinflag && strcmp(update->unit_style,"metal") != 0) + error->all(FLERR,"'metal' units have to be used with spins"); + + if (slabflag == 0 && domain->nonperiodic > 0) + error->all(FLERR,"Cannot use nonperiodic boundaries with PPPMDipoleSpin"); + if (slabflag) { + if (domain->xperiodic != 1 || domain->yperiodic != 1 || + domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1) + error->all(FLERR,"Incorrect boundaries with slab PPPMDipoleSpin"); + } + + if (order < 2 || order > MAXORDER) { + char str[128]; + sprintf(str,"PPPMDipoleSpin order cannot be < 2 or > than %d",MAXORDER); + error->all(FLERR,str); + } + + // extract short-range Coulombic cutoff from pair style + + triclinic = domain->triclinic; + if (triclinic) + error->all(FLERR,"Cannot yet use triclinic cells with PPPMDipoleSpin"); + + pair_check(); + + int itmp = 0; + double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp); + if (p_cutoff == NULL) + error->all(FLERR,"KSpace style is incompatible with Pair style"); + cutoff = *p_cutoff; + + // kspace TIP4P not yet supported + + if (tip4pflag) + error->all(FLERR,"Cannot yet use TIP4P with PPPMDipoleSpin"); + + scale = 1.0; + hbar = force->hplanck/MY_2PI; // eV/(rad.THz) + mub = 5.78901e-5; // in eV/T + mu_0 = 1.2566370614e-6; // in T.m/A + mub2mu0 = mub * mub * mu_0; // in eV + mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz + spsum_spsq(); + natoms_original = atom->natoms; + + // set accuracy (force units) from accuracy_relative or accuracy_absolute + + // is two_charge_force still relevant for spin systems? + + if (accuracy_absolute >= 0.0) accuracy = accuracy_absolute; + else accuracy = accuracy_relative * two_charge_force; + + // free all arrays previously allocated + + deallocate(); + if (peratom_allocate_flag) deallocate_peratom(); + + // setup FFT grid resolution and g_ewald + // normally one iteration thru while loop is all that is required + // if grid stencil does not extend beyond neighbor proc + // or overlap is allowed, then done + // else reduce order and try again + + int (*procneigh)[2] = comm->procneigh; + + GridComm *cgtmp = NULL; + int iteration = 0; + + while (order >= minorder) { + if (iteration && me == 0) + error->warning(FLERR,"Reducing PPPMDipoleSpin order b/c stencil extends " + "beyond nearest neighbor processor"); + + compute_gf_denom(); + set_grid_global(sp2); + set_grid_local(); + if (overlap_allowed) break; + + cgtmp = new GridComm(lmp,world,1,1, + nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, + nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out, + procneigh[0][0],procneigh[0][1],procneigh[1][0], + procneigh[1][1],procneigh[2][0],procneigh[2][1]); + cgtmp->ghost_notify(); + if (!cgtmp->ghost_overlap()) break; + delete cgtmp; + + order--; + iteration++; + } + + if (order < minorder) error->all(FLERR,"PPPMDipoleSpin order < minimum allowed order"); + if (!overlap_allowed && cgtmp->ghost_overlap()) + error->all(FLERR,"PPPMDipoleSpin grid stencil extends " + "beyond nearest neighbor processor"); + if (cgtmp) delete cgtmp; + + // adjust g_ewald + + if (!gewaldflag) adjust_gewald(); + + // calculate the final accuracy + + double estimated_accuracy = final_accuracy_dipole(sp2); + + // print stats + + int ngrid_max,nfft_both_max; + MPI_Allreduce(&ngrid,&ngrid_max,1,MPI_INT,MPI_MAX,world); + MPI_Allreduce(&nfft_both,&nfft_both_max,1,MPI_INT,MPI_MAX,world); + + if (me == 0) { + +#ifdef FFT_SINGLE + const char fft_prec[] = "single"; +#else + const char fft_prec[] = "double"; +#endif + + if (screen) { + fprintf(screen," G vector (1/distance) = %g\n",g_ewald); + fprintf(screen," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm); + fprintf(screen," stencil order = %d\n",order); + fprintf(screen," estimated absolute RMS force accuracy = %g\n", + estimated_accuracy); + fprintf(screen," estimated relative force accuracy = %g\n", + estimated_accuracy/two_charge_force); + fprintf(screen," using %s precision FFTs\n",fft_prec); + fprintf(screen," 3d grid and FFT values/proc = %d %d\n", + ngrid_max,nfft_both_max); + } + if (logfile) { + fprintf(logfile," G vector (1/distance) = %g\n",g_ewald); + fprintf(logfile," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm); + fprintf(logfile," stencil order = %d\n",order); + fprintf(logfile," estimated absolute RMS force accuracy = %g\n", + estimated_accuracy); + fprintf(logfile," estimated relative force accuracy = %g\n", + estimated_accuracy/two_charge_force); + fprintf(logfile," using %s precision FFTs\n",fft_prec); + fprintf(logfile," 3d grid and FFT values/proc = %d %d\n", + ngrid_max,nfft_both_max); + } + } + + // allocate K-space dependent memory + // don't invoke allocate peratom(), will be allocated when needed + + allocate(); + cg_dipole->ghost_notify(); + cg_dipole->setup(); + + // pre-compute Green's function denominator expansion + // pre-compute 1d charge distribution coefficients + + compute_gf_denom(); + compute_rho_coeff(); +} + +/* ---------------------------------------------------------------------- + compute the PPPMDipoleSpin long-range force, energy, virial +------------------------------------------------------------------------- */ + +void PPPMDipoleSpin::compute(int eflag, int vflag) +{ + int i,j; + + // set energy/virial flags + // invoke allocate_peratom() if needed for first time + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = evflag_atom = eflag_global = vflag_global = + eflag_atom = vflag_atom = 0; + + if (evflag_atom && !peratom_allocate_flag) { + allocate_peratom(); + cg_peratom_dipole->ghost_notify(); + cg_peratom_dipole->setup(); + } + + // if atom count has changed, update qsum and qsqsum + + if (atom->natoms != natoms_original) { + spsum_spsq(); + natoms_original = atom->natoms; + } + + // return if there are no spins + + if (spsqsum == 0.0) return; + + // convert atoms from box to lamda coords + + boxlo = domain->boxlo; + + // extend size of per-atom arrays if necessary + + if (atom->nmax > nmax) { + memory->destroy(part2grid); + nmax = atom->nmax; + memory->create(part2grid,nmax,3,"pppm_spin:part2grid"); + } + + // find grid points for all my particles + // map my particle charge onto my local 3d on-grid density + + particle_map(); + make_rho_spin(); + + // all procs communicate density values from their ghost cells + // to fully sum contribution in their 3d bricks + // remap from 3d decomposition to FFT decomposition + + cg_dipole->reverse_comm(this,REVERSE_SP); + brick2fft_dipole(); + + // compute potential gradient on my FFT grid and + // portion of e_long on this proc's FFT grid + // return gradients (electric fields) in 3d brick decomposition + // also performs per-atom calculations via poisson_peratom() + + poisson_ik_dipole(); + + // all procs communicate E-field values + // to fill ghost cells surrounding their 3d bricks + + cg_dipole->forward_comm(this,FORWARD_SP); + + // extra per-atom energy/virial communication + + if (evflag_atom) { + cg_peratom_dipole->forward_comm(this,FORWARD_SP_PERATOM); + } + + // calculate the force on my particles + + fieldforce_ik_spin(); + + // extra per-atom energy/virial communication + + if (evflag_atom) fieldforce_peratom_spin(); + + // sum global energy across procs and add in volume-dependent term + + const double spscale = mub2mu0 * scale; + const double g3 = g_ewald*g_ewald*g_ewald; + + if (eflag_global) { + double energy_all; + MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); + energy = energy_all; + + energy *= 0.5*volume; + energy -= spsqsum*2.0*g3/3.0/MY_PIS; + energy *= spscale; + } + + // sum global virial across procs + + if (vflag_global) { + double virial_all[6]; + MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world); + for (i = 0; i < 6; i++) virial[i] = 0.5*spscale*volume*virial_all[i]; + } + + // per-atom energy/virial + // energy includes self-energy correction + + if (evflag_atom) { + double **sp = atom->sp; + double spx,spy,spz; + int nlocal = atom->nlocal; + int ntotal = nlocal; + + if (eflag_atom) { + for (i = 0; i < nlocal; i++) { + spx = sp[i][0]*sp[i][3]; + spy = sp[i][1]*sp[i][3]; + spz = sp[i][2]*sp[i][3]; + eatom[i] *= 0.5; + eatom[i] -= (spx*spx + spy*spy + spz*spz)*2.0*g3/3.0/MY_PIS; + eatom[i] *= spscale; + } + } + + if (vflag_atom) { + for (i = 0; i < ntotal; i++) + for (j = 0; j < 6; j++) vatom[i][j] *= 0.5*spscale; + } + } + + // 2d slab correction + + if (slabflag == 1) slabcorr(); +} + +/* ---------------------------------------------------------------------- + create discretized "density" on section of global grid due to my particles + density(x,y,z) = charge "density" at grid points of my 3d brick + (nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts) + in global grid +------------------------------------------------------------------------- */ + +void PPPMDipoleSpin::make_rho_spin() +{ + int l,m,n,nx,ny,nz,mx,my,mz; + FFT_SCALAR dx,dy,dz; + FFT_SCALAR x0,y0,z0; + FFT_SCALAR x1,y1,z1; + FFT_SCALAR x2,y2,z2; + + // clear 3d density array + + memset(&(densityx_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0, + ngrid*sizeof(FFT_SCALAR)); + memset(&(densityy_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0, + ngrid*sizeof(FFT_SCALAR)); + memset(&(densityz_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0, + ngrid*sizeof(FFT_SCALAR)); + + // loop over my charges, add their contribution to nearby grid points + // (nx,ny,nz) = global coords of grid pt to "lower left" of charge + // (dx,dy,dz) = distance to "lower left" grid pt + // (mx,my,mz) = global coords of moving stencil pt + + double **sp = atom->sp; + double spx,spy,spz; + double **x = atom->x; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + + nx = part2grid[i][0]; + ny = part2grid[i][1]; + nz = part2grid[i][2]; + dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv; + dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv; + dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv; + + compute_rho1d(dx,dy,dz); + + spx = sp[i][0]*sp[i][3]; + spy = sp[i][1]*sp[i][3]; + spz = sp[i][2]*sp[i][3]; + z0 = delvolinv * spx; + z1 = delvolinv * spy; + z2 = delvolinv * spz; + for (n = nlower; n <= nupper; n++) { + mz = n+nz; + y0 = z0*rho1d[2][n]; + y1 = z1*rho1d[2][n]; + y2 = z2*rho1d[2][n]; + for (m = nlower; m <= nupper; m++) { + my = m+ny; + x0 = y0*rho1d[1][m]; + x1 = y1*rho1d[1][m]; + x2 = y2*rho1d[1][m]; + for (l = nlower; l <= nupper; l++) { + mx = l+nx; + densityx_brick_dipole[mz][my][mx] += x0*rho1d[0][l]; + densityy_brick_dipole[mz][my][mx] += x1*rho1d[0][l]; + densityz_brick_dipole[mz][my][mx] += x2*rho1d[0][l]; + } + } + } + } +} + +/* ---------------------------------------------------------------------- + interpolate from grid to get magnetic field & force on my particles for ik +------------------------------------------------------------------------- */ + +void PPPMDipoleSpin::fieldforce_ik_spin() +{ + int i,l,m,n,nx,ny,nz,mx,my,mz; + FFT_SCALAR dx,dy,dz; + FFT_SCALAR x0,y0,z0; + FFT_SCALAR ex,ey,ez; + FFT_SCALAR vxx,vyy,vzz,vxy,vxz,vyz; + + // loop over my charges, interpolate electric field from nearby grid points + // (nx,ny,nz) = global coords of grid pt to "lower left" of charge + // (dx,dy,dz) = distance to "lower left" grid pt + // (mx,my,mz) = global coords of moving stencil pt + + double **sp = atom->sp; + double spx,spy,spz; + double **x = atom->x; + double **f = atom->f; + double **fm = atom->fm; + + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) { + nx = part2grid[i][0]; + ny = part2grid[i][1]; + nz = part2grid[i][2]; + dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv; + dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv; + dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv; + + compute_rho1d(dx,dy,dz); + + ex = ey = ez = ZEROF; + vxx = vyy = vzz = vxy = vxz = vyz = ZEROF; + for (n = nlower; n <= nupper; n++) { + mz = n+nz; + z0 = rho1d[2][n]; + for (m = nlower; m <= nupper; m++) { + my = m+ny; + y0 = z0*rho1d[1][m]; + for (l = nlower; l <= nupper; l++) { + mx = l+nx; + x0 = y0*rho1d[0][l]; + ex -= x0*ux_brick_dipole[mz][my][mx]; + ey -= x0*uy_brick_dipole[mz][my][mx]; + ez -= x0*uz_brick_dipole[mz][my][mx]; + vxx -= x0*vdxx_brick_dipole[mz][my][mx]; + vyy -= x0*vdyy_brick_dipole[mz][my][mx]; + vzz -= x0*vdzz_brick_dipole[mz][my][mx]; + vxy -= x0*vdxy_brick_dipole[mz][my][mx]; + vxz -= x0*vdxz_brick_dipole[mz][my][mx]; + vyz -= x0*vdyz_brick_dipole[mz][my][mx]; + } + } + } + + // convert M-field to mech. and mag. forces + + const double spfactor = mub2mu0 * scale; + spx = sp[i][0]*sp[i][3]; + spy = sp[i][1]*sp[i][3]; + spz = sp[i][2]*sp[i][3]; + f[i][0] += spfactor*(vxx*spx + vxy*spy + vxz*spz); + f[i][1] += spfactor*(vxy*spx + vyy*spy + vyz*spz); + f[i][2] += spfactor*(vxz*spx + vyz*spy + vzz*spz); + + const double spfactorh = mub2mu0hbinv * scale; + fm[i][0] += spfactorh*ex; + fm[i][1] += spfactorh*ey; + fm[i][2] += spfactorh*ez; + + // create a new vector (in atom_spin style ?) to store long-range fm tables + + } +} + +/* ---------------------------------------------------------------------- + interpolate from grid to get per-atom energy/virial +------------------------------------------------------------------------- */ + +void PPPMDipoleSpin::fieldforce_peratom_spin() +{ + int i,l,m,n,nx,ny,nz,mx,my,mz; + FFT_SCALAR dx,dy,dz,x0,y0,z0; + FFT_SCALAR ux,uy,uz; + FFT_SCALAR v0x,v1x,v2x,v3x,v4x,v5x; + FFT_SCALAR v0y,v1y,v2y,v3y,v4y,v5y; + FFT_SCALAR v0z,v1z,v2z,v3z,v4z,v5z; + + // loop over my charges, interpolate from nearby grid points + // (nx,ny,nz) = global coords of grid pt to "lower left" of charge + // (dx,dy,dz) = distance to "lower left" grid pt + // (mx,my,mz) = global coords of moving stencil pt + + double **sp = atom->sp; + double spx,spy,spz; + double **x = atom->x; + + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) { + nx = part2grid[i][0]; + ny = part2grid[i][1]; + nz = part2grid[i][2]; + dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv; + dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv; + dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv; + + compute_rho1d(dx,dy,dz); + + ux = uy = uz = ZEROF; + v0x = v1x = v2x = v3x = v4x = v5x = ZEROF; + v0y = v1y = v2y = v3y = v4y = v5y = ZEROF; + v0z = v1z = v2z = v3z = v4z = v5z = ZEROF; + for (n = nlower; n <= nupper; n++) { + mz = n+nz; + z0 = rho1d[2][n]; + for (m = nlower; m <= nupper; m++) { + my = m+ny; + y0 = z0*rho1d[1][m]; + for (l = nlower; l <= nupper; l++) { + mx = l+nx; + x0 = y0*rho1d[0][l]; + if (eflag_atom) { + ux += x0*ux_brick_dipole[mz][my][mx]; + uy += x0*uy_brick_dipole[mz][my][mx]; + uz += x0*uz_brick_dipole[mz][my][mx]; + } + if (vflag_atom) { + v0x += x0*v0x_brick_dipole[mz][my][mx]; + v1x += x0*v1x_brick_dipole[mz][my][mx]; + v2x += x0*v2x_brick_dipole[mz][my][mx]; + v3x += x0*v3x_brick_dipole[mz][my][mx]; + v4x += x0*v4x_brick_dipole[mz][my][mx]; + v5x += x0*v5x_brick_dipole[mz][my][mx]; + v0y += x0*v0y_brick_dipole[mz][my][mx]; + v1y += x0*v1y_brick_dipole[mz][my][mx]; + v2y += x0*v2y_brick_dipole[mz][my][mx]; + v3y += x0*v3y_brick_dipole[mz][my][mx]; + v4y += x0*v4y_brick_dipole[mz][my][mx]; + v5y += x0*v5y_brick_dipole[mz][my][mx]; + v0z += x0*v0z_brick_dipole[mz][my][mx]; + v1z += x0*v1z_brick_dipole[mz][my][mx]; + v2z += x0*v2z_brick_dipole[mz][my][mx]; + v3z += x0*v3z_brick_dipole[mz][my][mx]; + v4z += x0*v4z_brick_dipole[mz][my][mx]; + v5z += x0*v5z_brick_dipole[mz][my][mx]; + } + } + } + } + + spx = sp[i][0]*sp[i][3]; + spy = sp[i][1]*sp[i][3]; + spz = sp[i][2]*sp[i][3]; + if (eflag_atom) eatom[i] += spx*ux + spy*uy + spz*uz; + if (vflag_atom) { + vatom[i][0] += spx*v0x + spy*v0y + spz*v0z; + vatom[i][1] += spx*v1x + spy*v1y + spz*v1z; + vatom[i][2] += spx*v2x + spy*v2y + spz*v2z; + vatom[i][3] += spx*v3x + spy*v3y + spz*v3z; + vatom[i][4] += spx*v4x + spy*v4y + spz*v4z; + vatom[i][5] += spx*v5x + spy*v5y + spz*v5z; + } + } +} + +/* ---------------------------------------------------------------------- + Slab-geometry correction term to dampen inter-slab interactions between + periodically repeating slabs. Yields good approximation to 2D Ewald if + adequate empty space is left between repeating slabs (J. Chem. Phys. + 111, 3155). Slabs defined here to be parallel to the xy plane. Also + extended to non-neutral systems (J. Chem. Phys. 131, 094107). +------------------------------------------------------------------------- */ + +void PPPMDipoleSpin::slabcorr() +{ + // compute local contribution to global spin moment + + double **x = atom->x; + double zprd = domain->zprd; + int nlocal = atom->nlocal; + + double spin = 0.0; + double **sp = atom->sp; + double spx,spy,spz; + for (int i = 0; i < nlocal; i++) { + spz = sp[i][2]*sp[i][3]; + spin += spz; + } + + // sum local contributions to get global spin moment + + double spin_all; + MPI_Allreduce(&spin,&spin_all,1,MPI_DOUBLE,MPI_SUM,world); + + // compute corrections + + const double e_slabcorr = MY_2PI*(spin_all*spin_all/12.0)/volume; + const double spscale = mub2mu0 * scale; + + if (eflag_global) energy += spscale * e_slabcorr; + + // per-atom energy + + if (eflag_atom) { + double efact = spscale * MY_2PI/volume/12.0; + for (int i = 0; i < nlocal; i++) { + spz = sp[i][2]*sp[i][3]; + eatom[i] += efact * spz * spin_all; + } + } + + // add on mag. force corrections + + double ffact = spscale * (-4.0*MY_PI/volume); + double **fm = atom->fm; + for (int i = 0; i < nlocal; i++) { + fm[i][2] += ffact * spin_all; + } +} + +/* ---------------------------------------------------------------------- + compute spsum,spsqsum,sp2 + called initially, when particle count changes, when spins are changed +------------------------------------------------------------------------- */ + +void PPPMDipoleSpin::spsum_spsq() +{ + const int nlocal = atom->nlocal; + + spsum = spsqsum = sp2 = 0.0; + if (atom->sp_flag) { + double **sp = atom->sp; + double spx, spy, spz; + double spsum_local(0.0), spsqsum_local(0.0); + + // not exactly the good loop: need to add norm of spins + + for (int i = 0; i < nlocal; i++) { + spx = sp[i][0]*sp[i][3]; + spy = sp[i][1]*sp[i][3]; + spz = sp[i][2]*sp[i][3]; + spsum_local += spx + spy + spz; + spsqsum_local += spx*spx + spy*spy + spz*spz; + } + + MPI_Allreduce(&spsum_local,&spsum,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&spsqsum_local,&spsqsum,1,MPI_DOUBLE,MPI_SUM,world); + + sp2 = spsqsum * mub2mu0; + } + + if (sp2 == 0 && comm->me == 0) + error->all(FLERR,"Using kspace solver PPPMDipoleSpin on system with no spins"); +} diff --git a/src/KSPACE/pppm_spin.h b/src/KSPACE/pppm_dipole_spin.h similarity index 66% rename from src/KSPACE/pppm_spin.h rename to src/KSPACE/pppm_dipole_spin.h index 3b4d42d4ea..347006e586 100644 --- a/src/KSPACE/pppm_spin.h +++ b/src/KSPACE/pppm_dipole_spin.h @@ -13,28 +13,23 @@ #ifdef KSPACE_CLASS -KSpaceStyle(pppm/spin,PPPMSpin) +KSpaceStyle(pppm/dipole/spin,PPPMDipoleSpin) #else -#ifndef LMP_PPPM_DIPOLE_H -#define LMP_PPPM_DIPOLE_H +#ifndef LMP_PPPM_DIPOLE_SPIN_H +#define LMP_PPPM_DIPOLE_SPIN_H -#include "pppm.h" +#include "pppm_dipole.h" namespace LAMMPS_NS { -class PPPMSpin : public PPPM { +class PPPMDipoleSpin : public PPPMDipole { public: - PPPMSpin(class LAMMPS *, int, char **); - virtual ~PPPMSpin(); + PPPMDipoleSpin(class LAMMPS *, int, char **); + virtual ~PPPMDipoleSpin(); void init(); - void setup(); - void setup_grid(); void compute(int, int); - int timing_1d(int, double &); - int timing_3d(int, double &); - double memory_usage(); protected: double hbar; // reduced Planck's constant @@ -42,55 +37,15 @@ class PPPMSpin : public PPPM { double mu_0; // vacuum permeability double mub2mu0; // prefactor for mech force double mub2mu0hbinv; // prefactor for mag force - void set_grid_global(); - double newton_raphson_f(); - - void allocate(); - void allocate_peratom(); - void deallocate(); - void deallocate_peratom(); - void compute_gf_denom(); void slabcorr(); - // grid communication - - void pack_forward(int, FFT_SCALAR *, int, int *); - void unpack_forward(int, FFT_SCALAR *, int, int *); - void pack_reverse(int, FFT_SCALAR *, int, int *); - void unpack_reverse(int, FFT_SCALAR *, int, int *); - // spin - FFT_SCALAR ***densityx_brick_spin,***densityy_brick_spin,***densityz_brick_spin; - FFT_SCALAR ***vdxx_brick_spin,***vdyy_brick_spin,***vdzz_brick_spin; - FFT_SCALAR ***vdxy_brick_spin,***vdxz_brick_spin,***vdyz_brick_spin; - FFT_SCALAR ***ux_brick_spin,***uy_brick_spin,***uz_brick_spin; - FFT_SCALAR ***v0x_brick_spin,***v1x_brick_spin,***v2x_brick_spin; - FFT_SCALAR ***v3x_brick_spin,***v4x_brick_spin,***v5x_brick_spin; - FFT_SCALAR ***v0y_brick_spin,***v1y_brick_spin,***v2y_brick_spin; - FFT_SCALAR ***v3y_brick_spin,***v4y_brick_spin,***v5y_brick_spin; - FFT_SCALAR ***v0z_brick_spin,***v1z_brick_spin,***v2z_brick_spin; - FFT_SCALAR ***v3z_brick_spin,***v4z_brick_spin,***v5z_brick_spin; - FFT_SCALAR *work3,*work4; - FFT_SCALAR *densityx_fft_spin,*densityy_fft_spin,*densityz_fft_spin; - class GridComm *cg_spin; - class GridComm *cg_peratom_spin; - int only_spin_flag; double spsum,spsqsum,sp2; - double find_gewald_spin(double, double, bigint, double, double); - double newton_raphson_f_spin(double, double, bigint, double, double); - double derivf_spin(double, double, bigint, double, double); - double compute_df_kspace_spin(); - double compute_qopt_spin(); - void compute_gf_spin(); void make_rho_spin(); - void brick2fft_spin(); - void poisson_ik_spin(); - void poisson_peratom_spin(); void fieldforce_ik_spin(); void fieldforce_peratom_spin(); - double final_accuracy_spin(); void spsum_spsq(); }; @@ -102,9 +57,9 @@ class PPPMSpin : public PPPM { /* ERROR/WARNING messages: -E: Cannot (yet) use charges with Kspace style PPPMSpin +E: Cannot (yet) use charges with Kspace style PPPMDipoleSpin -Charge-spin interactions are not yet implemented in PPPMSpin so this +Charge-spin interactions are not yet implemented in PPPMDipoleSpin so this feature is not yet supported. E: Must redefine kspace_style after changing to triclinic box @@ -123,11 +78,11 @@ E: Cannot (yet) use 'electron' units with spins This feature is not yet supported. -E: Cannot yet use triclinic cells with PPPMSpin +E: Cannot yet use triclinic cells with PPPMDipoleSpin This feature is not yet supported. -E: Cannot yet use TIP4P with PPPMSpin +E: Cannot yet use TIP4P with PPPMDipoleSpin This feature is not yet supported. @@ -207,9 +162,9 @@ outside a processor's sub-domain or even the entire simulation box. This indicates bad physics, e.g. due to highly overlapping atoms, too large a timestep, etc. -E: Using kspace solver PPPMSpin on system with no spins +E: Using kspace solver PPPMDipoleSpin on system with no spins -Must have non-zero spins with PPPMSpin. +Must have non-zero spins with PPPMDipoleSpin. E: Must use kspace_modify gewald for system with no spins diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp index 6460a6185f..fb2b6dd797 100644 --- a/src/SPIN/atom_vec_spin.cpp +++ b/src/SPIN/atom_vec_spin.cpp @@ -48,7 +48,7 @@ AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp) comm_x_only = 0; comm_f_only = 0; size_forward = 7; - size_reverse = 6; + size_reverse = 9; size_border = 10; size_velocity = 3; size_data_atom = 9; @@ -58,7 +58,6 @@ AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp) atom->sp_flag = 1; } - /* ---------------------------------------------------------------------- grow atom arrays n = 0 grows arrays by a chunk @@ -88,6 +87,7 @@ void AtomVecSpin::grow(int n) sp = memory->grow(atom->sp,nmax,4,"atom:sp"); fm = memory->grow(atom->fm,nmax*comm->nthreads,3,"atom:fm"); + fm_long = memory->grow(atom->fm_long,nmax*comm->nthreads,3,"atom:fm_long"); if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -103,7 +103,7 @@ void AtomVecSpin::grow_reset() tag = atom->tag; type = atom->type; mask = atom->mask; image = atom->image; x = atom->x; v = atom->v; f = atom->f; - sp = atom->sp; fm = atom->fm; + sp = atom->sp; fm = atom->fm; fm_long = atom->fm_long; } @@ -342,6 +342,9 @@ int AtomVecSpin::pack_reverse(int n, int first, double *buf) buf[m++] = fm[i][0]; buf[m++] = fm[i][1]; buf[m++] = fm[i][2]; + buf[m++] = fm_long[i][0]; + buf[m++] = fm_long[i][1]; + buf[m++] = fm_long[i][2]; } return m; @@ -361,6 +364,9 @@ void AtomVecSpin::unpack_reverse(int n, int *list, double *buf) fm[j][0] += buf[m++]; fm[j][1] += buf[m++]; fm[j][2] += buf[m++]; + fm_long[j][0] += buf[m++]; + fm_long[j][1] += buf[m++]; + fm_long[j][2] += buf[m++]; } } @@ -939,6 +945,7 @@ bigint AtomVecSpin::memory_usage() if (atom->memcheck("sp")) bytes += memory->usage(sp,nmax,4); if (atom->memcheck("fm")) bytes += memory->usage(fm,nmax*comm->nthreads,3); + if (atom->memcheck("fm_long")) bytes += memory->usage(fm_long,nmax*comm->nthreads,3); return bytes; } @@ -947,6 +954,7 @@ void AtomVecSpin::force_clear(int n, size_t nbytes) { memset(&atom->f[0][0],0,3*nbytes); memset(&atom->fm[0][0],0,3*nbytes); + memset(&atom->fm_long[0][0],0,3*nbytes); } diff --git a/src/SPIN/atom_vec_spin.h b/src/SPIN/atom_vec_spin.h index 34bc55b99f..1f1c34df52 100644 --- a/src/SPIN/atom_vec_spin.h +++ b/src/SPIN/atom_vec_spin.h @@ -68,10 +68,12 @@ class AtomVecSpin : public AtomVec { int *type,*mask; imageint *image; double **x,**v,**f; // lattice quantities - double **sp,**fm; // spin quantities - // sp[i][0-2] direction of the spin i + + // spin quantities + double **sp; // sp[i][0-2] direction of the spin i // sp[i][3] atomic magnetic moment of the spin i - + double **fm; // fm[i][0-2] direction of magnetic precession + double **fm_long; // storage of long-range spin prec. components }; } diff --git a/src/SPIN/pair_spin_exchange.cpp b/src/SPIN/pair_spin_exchange.cpp index cc074bb97d..74570afbce 100644 --- a/src/SPIN/pair_spin_exchange.cpp +++ b/src/SPIN/pair_spin_exchange.cpp @@ -441,8 +441,6 @@ void PairSpinExchange::allocate() memory->create(cutsq,n+1,n+1,"pair:cutsq"); } - - /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ @@ -527,4 +525,3 @@ void PairSpinExchange::read_restart_settings(FILE *fp) MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } - diff --git a/src/SPIN/pair_spin_long.cpp b/src/SPIN/pair_spin_long.cpp index 95c4e6b5a9..ca4f3d5e24 100644 --- a/src/SPIN/pair_spin_long.cpp +++ b/src/SPIN/pair_spin_long.cpp @@ -13,19 +13,19 @@ /* ------------------------------------------------------------------------ Contributing authors: Julien Tranchida (SNL) - Stan Moore (SNL) - - Please cite the related publication: - Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). - Massively parallel symplectic algorithm for coupled magnetic spin dynamics - and molecular dynamics. arXiv preprint arXiv:1801.10233. -------------------------------------------------------------------------- */ + Aidan Thompson (SNL) + Please cite the related publication: + Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). + Massively parallel symplectic algorithm for coupled magnetic spin dynamics + and molecular dynamics. Journal of Computational Physics. +------------------------------------------------------------------------- */ #include #include #include #include + #include "pair_spin_long.h" #include "atom.h" #include "comm.h" @@ -80,10 +80,154 @@ PairSpinLong::~PairSpinLong() { if (allocated) { memory->destroy(setflag); - memory->destroy(cutsq); + memory->destroy(cut_spin_long); } } +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairSpinLong::settings(int narg, char **arg) +{ + if (narg < 1 || narg > 2) + error->all(FLERR,"Incorrect args in pair_style command"); + + if (strcmp(update->unit_style,"metal") != 0) + error->all(FLERR,"Spin simulations require metal unit style"); + + cut_spin_long_global = force->numeric(FLERR,arg[0]); + //cut_spin = force->numeric(FLERR,arg[0]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) { + for (j = i+1; j <= atom->ntypes; j++) { + if (setflag[i][j]) { + cut_spin_long[i][j] = cut_spin_long_global; + } + } + } + } + +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairSpinLong::coeff(int narg, char **arg) +{ + if (!allocated) allocate(); + + // check if args correct + + if (strcmp(arg[2],"long") != 0) + error->all(FLERR,"Incorrect args in pair_style command"); + if (narg != 3) + error->all(FLERR,"Incorrect args in pair_style command"); + + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + double spin_long_cut_one = force->numeric(FLERR,arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + setflag[i][j] = 1; + cut_spin_long[i][j] = spin_long_cut_one; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairSpinLong::init_style() +{ + if (!atom->sp_flag) + error->all(FLERR,"Pair spin requires atom/spin style"); + + // need a full neighbor list + + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + + // checking if nve/spin is a listed fix + + int ifix = 0; + while (ifix < modify->nfix) { + if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break; + ifix++; + } + if (ifix == modify->nfix) + error->all(FLERR,"pair/spin style requires nve/spin"); + + // get the lattice_flag from nve/spin + + for (int i = 0; i < modify->nfix; i++) { + if (strcmp(modify->fix[i]->style,"nve/spin") == 0) { + lockfixnvespin = (FixNVESpin *) modify->fix[i]; + lattice_flag = lockfixnvespin->lattice_flag; + } + } + + // insure use of KSpace long-range solver, set g_ewald + + if (force->kspace == NULL) + error->all(FLERR,"Pair style requires a KSpace style"); + + g_ewald = force->kspace->g_ewald; + +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairSpinLong::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + + cut_spin_long[j][i] = cut_spin_long[i][j]; + + return cut_spin_long_global; +} + +/* ---------------------------------------------------------------------- + extract the larger cutoff if "cut" or "cut_coul" +------------------------------------------------------------------------- */ + +void *PairSpinLong::extract(const char *str, int &dim) +{ + if (strcmp(str,"cut") == 0) { + dim = 0; + return (void *) &cut_spin_long_global; + } else if (strcmp(str,"cut_coul") == 0) { + dim = 0; + return (void *) &cut_spin_long_global; + } else if (strcmp(str,"ewald_order") == 0) { + ewald_order = 0; + ewald_order |= 1<<1; + ewald_order |= 1<<3; + dim = 0; + return (void *) &ewald_order; + } else if (strcmp(str,"ewald_mix") == 0) { + dim = 0; + return (void *) &mix_flag; + } + return NULL; +} + /* ---------------------------------------------------------------------- */ void PairSpinLong::compute(int eflag, int vflag) @@ -95,6 +239,7 @@ void PairSpinLong::compute(int eflag, int vflag) double evdwl,ecoul; double xi[3],rij[3]; double spi[4],spj[4],fi[3],fmi[3]; + double local_cut2; double pre1,pre2,pre3; int *ilist,*jlist,*numneigh,**firstneigh; @@ -156,26 +301,25 @@ void PairSpinLong::compute(int eflag, int vflag) rij[2] = x[j][2] - xi[2]; rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]; - if (rsq < cutsq[itype][jtype]) { + local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype]; + + if (rsq < local_cut2) { r2inv = 1.0/rsq; rinv = sqrt(r2inv); - if (rsq < cut_spinsq) { - r = sqrt(rsq); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - bij[0] = erfc * rinv; - bij[1] = (bij[0] + pre1*expm2) * r2inv; - bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv; - bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv; + bij[0] = erfc * rinv; + bij[1] = (bij[0] + pre1*expm2) * r2inv; + bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv; + bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv; - compute_long(i,j,rij,bij,fmi,spi,spj); - compute_long_mech(i,j,rij,bij,fmi,spi,spj); - - } + compute_long(i,j,rij,bij,fmi,spi,spj); + compute_long_mech(i,j,rij,bij,fmi,spi,spj); } // force accumulation @@ -194,7 +338,7 @@ void PairSpinLong::compute(int eflag, int vflag) } if (eflag) { - if (rsq <= cut_spinsq) { + if (rsq <= local_cut2) { evdwl -= spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]; evdwl *= hbar; @@ -219,6 +363,7 @@ void PairSpinLong::compute_single_pair(int ii, double fmi[3]) double r,rinv,r2inv,rsq; double grij,expm2,t,erfc; double bij[4],xi[3],rij[3],spi[4],spj[4]; + double local_cut2; double pre1,pre2,pre3; int *ilist,*jlist,*numneigh,**firstneigh; @@ -267,25 +412,24 @@ void PairSpinLong::compute_single_pair(int ii, double fmi[3]) rij[2] = x[j][2] - xi[2]; rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]; - if (rsq < cutsq[itype][jtype]) { + local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype]; + + if (rsq < local_cut2) { r2inv = 1.0/rsq; rinv = sqrt(r2inv); - if (rsq < cut_spinsq) { - r = sqrt(rsq); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - bij[0] = erfc * rinv; - bij[1] = (bij[0] + pre1*expm2) * r2inv; - bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv; - bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv; + bij[0] = erfc * rinv; + bij[1] = (bij[0] + pre1*expm2) * r2inv; + bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv; + bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv; - compute_long(i,j,rij,bij,fmi,spi,spj); - - } + compute_long(i,j,rij,bij,fmi,spi,spj); } } @@ -361,111 +505,7 @@ void PairSpinLong::allocate() for (int j = i; j <= n; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairSpinLong::settings(int narg, char **arg) -{ - if (narg < 1 || narg > 2) - error->all(FLERR,"Incorrect args in pair_style command"); - - if (strcmp(update->unit_style,"metal") != 0) - error->all(FLERR,"Spin simulations require metal unit style"); - - cut_spin = force->numeric(FLERR,arg[0]); - -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairSpinLong::coeff(int narg, char **arg) -{ - if (narg < 4 || narg > 5) - error->all(FLERR,"Incorrect args for pair coefficients"); - if (!allocated) allocate(); - - // check if args correct - - if (strcmp(arg[2],"long") != 0) - error->all(FLERR,"Incorrect args in pair_style command"); - if (narg != 3) - error->all(FLERR,"Incorrect args in pair_style command"); - - int ilo,ihi,jlo,jhi; - force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); - force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - setflag[i][j] = 1; - count++; - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairSpinLong::init_one(int i, int j) -{ - if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); - - double cut = cut_spin; - return cut; -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairSpinLong::init_style() -{ - if (!atom->sp_flag) - error->all(FLERR,"Pair spin requires atom/spin style"); - - // need a full neighbor list - - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; - - // checking if nve/spin is a listed fix - - int ifix = 0; - while (ifix < modify->nfix) { - if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break; - ifix++; - } - if (ifix == modify->nfix) - error->all(FLERR,"pair/spin style requires nve/spin"); - - // get the lattice_flag from nve/spin - - for (int i = 0; i < modify->nfix; i++) { - if (strcmp(modify->fix[i]->style,"nve/spin") == 0) { - lockfixnvespin = (FixNVESpin *) modify->fix[i]; - lattice_flag = lockfixnvespin->lattice_flag; - } - } - - // insure use of KSpace long-range solver, set g_ewald - - if (force->kspace == NULL) - error->all(FLERR,"Pair style requires a KSpace style"); - - g_ewald = force->kspace->g_ewald; - - cut_spinsq = cut_spin * cut_spin; + memory->create(cut_spin_long,n+1,n+1,"pair:cut_spin_long"); } /* ---------------------------------------------------------------------- @@ -477,10 +517,14 @@ void PairSpinLong::write_restart(FILE *fp) write_restart_settings(fp); int i,j; - for (i = 1; i <= atom->ntypes; i++) + for (i = 1; i <= atom->ntypes; i++) { for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&cut_spin_long[i][j],sizeof(int),1,fp); + } } + } } /* ---------------------------------------------------------------------- @@ -495,11 +539,18 @@ void PairSpinLong::read_restart(FILE *fp) int i,j; int me = comm->me; - for (i = 1; i <= atom->ntypes; i++) + for (i = 1; i <= atom->ntypes; i++) { for (j = i; j <= atom->ntypes; j++) { if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&cut_spin_long[i][j],sizeof(int),1,fp); + } + MPI_Bcast(&cut_spin_long[i][j],1,MPI_INT,0,world); + } } + } } /* ---------------------------------------------------------------------- @@ -508,7 +559,7 @@ void PairSpinLong::read_restart(FILE *fp) void PairSpinLong::write_restart_settings(FILE *fp) { - fwrite(&cut_spin,sizeof(double),1,fp); + fwrite(&cut_spin_long_global,sizeof(double),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } @@ -519,32 +570,9 @@ void PairSpinLong::write_restart_settings(FILE *fp) void PairSpinLong::read_restart_settings(FILE *fp) { if (comm->me == 0) { - fread(&cut_spin,sizeof(double),1,fp); + fread(&cut_spin_long_global,sizeof(double),1,fp); fread(&mix_flag,sizeof(int),1,fp); } - MPI_Bcast(&cut_spin,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_spin_long_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } - -/* ---------------------------------------------------------------------- */ - -void *PairSpinLong::extract(const char *str, int &dim) -{ - if (strcmp(str,"cut") == 0) { - dim = 0; - return (void *) &cut_spin; - } else if (strcmp(str,"cut_coul") == 0) { - dim = 0; - return (void *) &cut_spin; - } else if (strcmp(str,"ewald_order") == 0) { - ewald_order = 0; - ewald_order |= 1<<1; - ewald_order |= 1<<3; - dim = 0; - return (void *) &ewald_order; - } else if (strcmp(str,"ewald_mix") == 0) { - dim = 0; - return (void *) &mix_flag; - } - return NULL; -} diff --git a/src/SPIN/pair_spin_long.h b/src/SPIN/pair_spin_long.h index 867b771f74..0cdf6d2b80 100644 --- a/src/SPIN/pair_spin_long.h +++ b/src/SPIN/pair_spin_long.h @@ -49,6 +49,8 @@ class PairSpinLong : public PairSpin { void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); + + double cut_spin_long_global; // global long cutoff distance protected: double hbar; // reduced Planck's constant @@ -56,7 +58,8 @@ class PairSpinLong : public PairSpin { double mu_0; // vacuum permeability double mub2mu0; // prefactor for mech force double mub2mu0hbinv; // prefactor for mag force - double cut_spin, cut_spinsq; + + double **cut_spin_long; // cutoff distance long double g_ewald; int ewald_order;