diff -Naur lammps-28Sep09/doc/JPG/atc_nanotube.jpg lammps-29Sep09/doc/JPG/atc_nanotube.jpg --- lammps-28Sep09/doc/JPG/atc_nanotube.jpg 1969-12-31 17:00:00.000000000 -0700 +++ lammps-29Sep09/doc/JPG/atc_nanotube.jpg 2009-09-28 13:50:12.000000000 -0600 @@ -0,0 +1,56 @@ +JFIFHHExifMM*JR(iZHHcC  +   $.' ",#(7),01444'9=82<.342C  2!!222222222222222222222222222222222222222222222222221"M !1AQa"2q#$Brt356RSs%7Cu4DbcdT3!1AQ"aq#2B$3 ?2d%e*q#̚2s<zީeK1De=kFw +cARTOrNblOlƷ:kMj@1=6sҙ#5,u50qE>:ɍ_)@kPح~?-6%K])hUƮ!qj+U7fhdx|5ejvsl=J`k~r$蓚Ē̆‡^װ􅥍a6 iؙ!iu+*ZI +ug>UlITrH*H# SR0l?*9e#]lŠy)eD $dg8K}u 﫽 X4ҡRJyJ)UKb")JR"")JR"")JR"")JR"")JR"")JR"")JR""MؿZ?MQ >X~LLF 0FV@^w$7Y^[b~YtgPX̂~f4h6irSa ('|T~(ExfX7*"auzTLJ2@~%usVݗT[Ӣ]1qƌĜl2=(28sW +.݉Xn,#)|q<ܛ|Ru#;Gʙ )PRbwis咮$v;Vɇq]Q.4ƒk XE%u!crht3=qWU9 #AF9'6HޚJRo G[?d/ohimO.SȧJPQ뿍Wl|Vnj%ҧd,@'+;ª%@Os-R.;0254V?bR_)JR)JR)JR)JR)JR)JR)JR#۵k;3%! uX8sq^å &xij ~GN?>%2LvR7[ `2LA'_mtQ*gB Q';jagqėI8l .< qHlmxS6ﭚč*]mS“fۙ 5 +A +I6r ^ps$7 |/%'EWŕ@x\'lf5Uwdi3QU(G-uikHȑ%`ʲ $sTЯ#W(4H1N`ֿ[VJ_o(VZq>Ӕ,-A*tحIIOk3}Hh&5瞟:H!QYSɩ $*]u}Lqp:AIo ++^6i[52h4D G*&>UAMpۭ&E˜CaT+4)JR)JR)JR)JR)JR)JR)JR5eJ"qכhGi6흉 +:G+)'E[0mG*nu+VB[ȹ*jaTGil'x_wmxΆ>{C?>$n"]퉗ƋIgmԦ[W^`qmFmÍ-+A)' qI/?^>"5UvkUk@#wz֩tR?7P]MkV8Cⱑ2\{:oK&91zC ,c2qj'ZUk%^\OTC[uwLk1HfLAՒ'ajߪ,y1Xw-8z.(ΙTj=fi\W5hS4")_!@¸.$NOJ^'+* +B7"")JR"")JR"")JR"")JR"")JR""M6B}$ SY?W]μfI“Z?Iʰ:»v:&9vPq/Zr2=k.Ļ< #c BeU49#bCrWo$CHn2?x"okWBk?##늎nG&%ȭi*DHW!i9 +muh +HVqd  _ ŅXͅqi IFЕ6yZ掵k|y\-|QR-PZf+l)ĩNx Í"GQߺ Po"[CW&˘5ڄ +45rWķjiږMt];&. X6B/T!| {iJ`Q +樴ܔp1w}>;+d3bVʀ(#V#ܢȌ򿤢)4>@q>Bj +L 廇zE]mi+* T{6Z)Ķ9n 딟yWGmҙ\hehuGcӈlvzU-\n`8#8L|]-xyimsJZl07NՊJέߊ^,'268z uT6F_+?hU5DKdR9랔7`{ۧi()%F$7͏[E5KU,N&le]M^եMu|q\%-s +$d$wD!"BFR$1ͽa8mq$y$uMy|E*tj^w<^u'luٯuׄgdx~ꤟo;/2x) 9ZUs֥ej\JHʖn=k"L)Rq +I>_ZEDZ䐷!r9_N Qkikhnέ6/w?Z[40p҂Bz"qOOV)\<dxcjJj;q(p _ RK-<ҥvef7EVYnӛ祰s(RR{LhK+XVG1 yb1x1%$,rN2rqZ v6TU)̶v`6Jχ6$3)E!@VҭNĆ m;JUJRCM!#*R$:(`Y'$jj p<ҠMպ۬GN\H9QOVg-ltbt| +}Xi.q}zV*"0JT[*j1RD%(emzx\qdgU±'*l$*8oҸjdRE)JDO?>Tuzijkr%0k+R]I߱w3QvmYb+òT<ʾU2BZPeɛf%{B}BBBIiT!9LjQخMYO +cmpljҩ6'=3mVT[pthØ%h\ Vkսz@r˼ ޣ8sE$%6Z_[$n c95hްR͙v:~XJ?|e?}E담i3R`8!Hl[)$18+dBۮɉ*<¹k \~21բWܵFd/M̼"S<@NFqGQLIhvl +q J8P= o]Q"T}˃SO,૮8ڧ 𖝉5qDIa*,OJЦqo}\A1oǞ]ZB02I\5.R#H: P#cf1ySC~>+it +d,W\Qi,ϊT ;+އvt&5Dwsy>. HjӀ~£[3m4ޓE륳%8"BDfc\_鶄u̷SOxoҮbԉp[j'Uܥ :x6P)Z }R) `#)=|jT}~&/U`@cﭲBC QJRLzVzn3eY. .Ԃr3^?:ܬss$}p!<@ ۿ¨(a|AQ ss$A-$*9j${Ǖt& N.ar8GQMOut-zJ-+lua K2ӍڐŜӎՔ顦;d1Lo3%:D}cc~wWng[quHWA]% IO +dI{Lq>gdcd~07i'm5ń{á< +"58G<pL\{=2\gwVv-CfRRCqČ`x)blw$! +K{)c;dt'Cs\޵c%^}g_n Itw^&rgQJ ҀR;>fS~;uhPROȊy @J] +ў?֬-7c^;3c.6H8qA*͈ j57\MJ:ש1_-HRE].S^ beOhLGC?^HZVfZޙk[ͥt' ]R{V՚ +촲;ʙf6ENHmp9ejly);}g2`[nG/R`2sYm}rNM!dy@:HvB͇.}d)M :2 ֥?&OVۯcE\."|¦=bK=<,຀J}xҖWbII!-,b;c#QVe  ` +KNqH5k %.us2J;DbLfςAYmR':R꠾8pHB^,"<9&v9a?Mx4ƥE\=qa)]BО/}AKp] _f4wd*fkjaDq`Bk8$)$zjZSMLɎťpUEѠ(xb>UcޭmHqj3+!J )H#aa$2؀"老qJ/WKt;tf2mn5^UeͷT{_Њ;ʝ5Bi8⒔$eJQ֗P#yn/]릻k ̑";/ٗ@C)Q9:- u,Iu{`"hoڙBKxl!228)q[ +]p[y%i Tb`xjnVv#>-8?i85 ވeL ,!8ޫVZ\+ HRC{gIWZ3Mo8:mWh~=,1)ԇ +U!c9ڻt]K[*5qW[c1'@o|`{YoԒCzd-9 B)m ,}ё~U2~[sD]e32ue{.;ٝJ<+O+NZ^r#a @^x {s/3/GVQ78>8EQ)߬5V)ø䧌}FO.BoqI:W^3.XXxr~鮇%LF/{?$wWsH2!*oH#m^ǒ0(_k$4 9*Lm< MiccҰږ!CmK!p7_ ޽J5\ H<)7Vj +lH ܅DK%La<ʸ%|]zVҺݫn<K!-95_ш5҆JpI*8+J#fIRC{x|ڭ4P^_\+N%*ۂMμ7oԞv2,Hw+`}< d ++"2hYy;WF>CWd~{c#2VprwV߇į-fٗ+m2,7kSejv#G +TʂԷj%dppHN oyՉ8o6j=@en6iQPO1 Lc+`j5K3PӃURq"kjm +q0 +FO}U.|8P: +ֵиK[< %` x T8l^uդ m^FzËPqiuq5Enshur1 iK&u!sm8¤:GxT& +ͨK^,CykJ敚h:; Şi[l QI 'a^yޥ$^8c6U$t+<^@PBBkn5PN"Ga5shLKiBCb! m>%ݥ`>6m)=C}{WdPg-JÉC/6‚Uw؊cmIZ皸hUSI[q.^~jԓ-}B 8穀\RҤg֪`g |FI  FHxדQj{fJ%@GP0Rv +Iת\1ر@v4*A71@|+\+Pl:-vdH n{@J׉Z%Ϳ܇FomOx48ċe疒㣳Kd:Ê 0h`VzHv&u}5Ӧ֏ԙ梸RXZHy4J”@uQ֜4*)! AY +s.W W-Ҕ)H~;2[-Aꕤ(}JIVҊ%[t~~qG%͵$!U*E^]Z\+7B[d\@jGi-s9Ejwr}MS-6ؿEB]y_j;HFڼ߃i^-¨SY?_.$j\8Rר#nq$z\Rs6WЖXz3$!*gmYYsycKu%&US0!D$<.6|Rk&MApHBYhs)Dy +{$4+4bIF[.L1r.ABTwQ!!D(c =>U'onEmCEC# p00 ! A@Rvv+\jdYioX]r#j)r3ݹ;wxU2Y@7C IHpW3tEvΩkDR5*݉)=yxW9@a9#v++M28m)\iqE+J-.40<+[)s fZC@䝼V)~r%}ªAjq+T4+nrf)JE)JDul0R>/[H.tQT+"ȭT[);H)>V6lȴ5nulYB*q<{9>BiVv =^=|EN[O.J=>CTXܙjQa+JR)JR)\W4O3cHƒUaR,1YU S'y%A'nm3+*b:uY ]Չw+"r#wbQ Y!^_^ČYdwb}̝p ~ t2Kaj7ZIUV8Q%([%$kƋ*wT!\g^LS c!(\ұQJR)JR)JR)JR)JR)JR)JRsJD,C} E]7F.4dVWB)9I;5??WkgzjE15=qmR O\ +R,:Ô|ku$+cwsCÒRc6l`9UjLEҷ2͵̊Ơ8R@Ě Q1t >ӜRRr9- JЍ~.eP뇈Pȓ\CԆEVָkˌŎGTktdx:ҟ2I;,g5ػyp$uiI@N'$/fk[LBS$B{ #aX^|*,ɺAq5 bHR]xU,[^,S#$u~>$5}>i^=읿 Nzoj w z1Ռ:zJHse5еDRm.A?jwsP56a1 FNJAKy*`[_l +ZkTRdo֩T@lrz˙ؐ5WsT攥")JR"")JR"")JR"")JR"")JR"")JR"")JR' +kWiI)x6@1vT>dV=*kE&/!+HJд- A4]Vb}zPX#tͱ$Hsv#L.sU{<auը M2 +ЎZ9 +Z"ʓ [J0-lz{Rbbyl?BA+iFEʾyk놲\񕿀nw_pOH}.->UT\%:҄!BNG>&llZ92*gFei1Υ jIϺRGόґ~1>2-55 ӃUQ/0Ӄ-@>{:nt6$tzVٓ\7%A2%%H{|+'0MGvHГqn/` TuWlM%͡*hF;M5%ؕbC| FB1CU0$ŎF[R [ҔD攥")JR"")JR"")JR"")JR"")JR"")JR"")JR ދZQ(%"TI8sT 8[[Y?hi9ERiݺ<)lo#57Hs.P6jS)ӎCG=(h$=665] +|I6`6Jdߢm0uo;<؏'%^zVIK][v)<,,yҵ$JI.OY[+Sk_@II`MgeYb`l.Nd}/M% Z5<12!&VhpaIR@#$}+5p6{P/7~!X'&;O!JBr +Tp|:J\Vra".KU?j{cl|!28πmy+h2Z[U޵ԵGS Ă';U Lk̛L_59# Vfk3Description:

This fix creates a coupled finite element (FE) and molecular dynamics -(MD) simulation and/or an on-the-fly estimation of continuum fields. -Interscale operators are defined that construct continuum fields from +(MD) simulation and/or an on-the-fly estimation of continuum fields, +where a FE mesh is specified and overlaps the particles, something +like this: +

+
+
+

Interscale operators are defined that construct continuum fields from atomic data. Coupled simulations use FE projection approximated on a discrete field. Currently, coupling is restricted to thermal physics. The Hardy module can use either FE projection or integration Kernels @@ -47,12 +52,14 @@ thermostat can apply heat sources to the atoms that varies in space on the same scale as the FE element size. Meshes are not created automatically and must be specified on LAMMPS regions with prescribed -element sizes. Note that mesh computations and storage run in serial -(not parallelized) so performance will degrade when large element -counts are used. +element sizes.

-

Note that coupling and post-processing can be combined in the same -simulations using separate fix atc commands. +

Coupling and post-processing can be combined in the same simulations +using separate fix atc commands. +

+

Note that mesh computations and storage run in serial (not +parallelized) so performance will degrade when large element counts +are used.

For detailed exposition of the theory and algorithms implemented in this fix, please see the papers here and here. diff -Naur lammps-28Sep09/doc/fix_atc.txt lammps-29Sep09/doc/fix_atc.txt --- lammps-28Sep09/doc/fix_atc.txt 2009-09-24 12:06:18.000000000 -0600 +++ lammps-29Sep09/doc/fix_atc.txt 2009-09-28 13:50:12.000000000 -0600 @@ -28,7 +28,12 @@ [Description:] This fix creates a coupled finite element (FE) and molecular dynamics -(MD) simulation and/or an on-the-fly estimation of continuum fields. +(MD) simulation and/or an on-the-fly estimation of continuum fields, +where a FE mesh is specified and overlaps the particles, something +like this: + +:c,image(JPG/atc_nanotube.jpg) + Interscale operators are defined that construct continuum fields from atomic data. Coupled simulations use FE projection approximated on a discrete field. Currently, coupling is restricted to thermal physics. @@ -40,12 +45,14 @@ thermostat can apply heat sources to the atoms that varies in space on the same scale as the FE element size. Meshes are not created automatically and must be specified on LAMMPS regions with prescribed -element sizes. Note that mesh computations and storage run in serial -(not parallelized) so performance will degrade when large element -counts are used. +element sizes. + +Coupling and post-processing can be combined in the same simulations +using separate fix atc commands. -Note that coupling and post-processing can be combined in the same -simulations using separate fix atc commands. +Note that mesh computations and storage run in serial (not +parallelized) so performance will degrade when large element counts +are used. For detailed exposition of the theory and algorithms implemented in this fix, please see the papers "here"_#Wagner and "here"_#Zimmerman. diff -Naur lammps-28Sep09/doc/mass.html lammps-29Sep09/doc/mass.html --- lammps-28Sep09/doc/mass.html 2009-06-26 12:22:33.000000000 -0600 +++ lammps-29Sep09/doc/mass.html 2009-09-29 08:10:23.000000000 -0600 @@ -65,9 +65,9 @@

If you define a hybrid atom style which includes one (or more) sub-styles which require per-type mass and one (or more) -sub-styles which require per-atom mass, then you must define both. In -this case the per-type mass will be ignored; only the per-atom mass -will be used by LAMMPS. +sub-styles which require per-atom mass, then you must define both. +However, in this case the per-type mass will be ignored; only the +per-atom mass will be used by LAMMPS.

Restrictions:

diff -Naur lammps-28Sep09/doc/mass.txt lammps-29Sep09/doc/mass.txt --- lammps-28Sep09/doc/mass.txt 2009-06-26 12:22:33.000000000 -0600 +++ lammps-29Sep09/doc/mass.txt 2009-09-29 08:10:23.000000000 -0600 @@ -62,9 +62,9 @@ If you define a "hybrid atom style"_atom_style.html which includes one (or more) sub-styles which require per-type mass and one (or more) -sub-styles which require per-atom mass, then you must define both. In -this case the per-type mass will be ignored; only the per-atom mass -will be used by LAMMPS. +sub-styles which require per-atom mass, then you must define both. +However, in this case the per-type mass will be ignored; only the +per-atom mass will be used by LAMMPS. [Restrictions:] diff -Naur lammps-28Sep09/doc/min_modify.html lammps-29Sep09/doc/min_modify.html --- lammps-28Sep09/doc/min_modify.html 2009-02-25 17:18:18.000000000 -0700 +++ lammps-29Sep09/doc/min_modify.html 2009-09-29 08:10:23.000000000 -0600 @@ -32,30 +32,29 @@

Description:

This command sets parameters that affect the energy minimization -algorithms. The various settings may affect the convergence rate and -overall number of force evaluations required by a minimization, so -users can experiment with these parameters to tune their -minimizations. -

-

The minimization algorithms have an outer iteration (conjugate -gradient or steepest descent) and an inner iteration which is steps -along a one-dimensional line search in a particular search direction. -The dmax parameter is how far any atom can move in a single line -search in any dimension (x, y, or z). Thus a value of 0.1 in real -units means no atom will move further than 0.1 Angstroms -in a single outer iteration. This prevents highly overlapped atoms -from being moved long distances (e.g. through another atom) due to -large forces. -

-

The choice of line search algorithm can be selected via the line -keyword. The default backtracking search is very robust and should -always find a local energy minimum. However, it will "converge" when -it can no longer reduce the energy of the system. Individual atom -forces may still be larger than desired at this point, because the -energy change is measured as the difference of two large values -(energy before and energy after) and that difference may be smaller -than machine epsilon even if atoms could move in the gradient -direction to reduce forces further. +algorithms selected by the min_style command. The +various settings may affect the convergence rate and overall number of +force evaluations required by a minimization, so users can experiment +with these parameters to tune their minimizations. +

+

The cg and sd minimization styles have an outer iteration and an +inner iteration which is steps along a one-dimensional line search in +a particular search direction. The dmax parameter is how far any +atom can move in a single line search in any dimension (x, y, or z). +Thus a value of 0.1 in real units means no atom will move +further than 0.1 Angstroms in a single outer iteration. This prevents +highly overlapped atoms from being moved long distances (e.g. through +another atom) due to large forces. +

+

The choice of line search algorithm for the cg and sd minimization +styles can be selected via the line keyword. The default +backtracking search is robust and should always find a local energy +minimum. However, it will "converge" when it can no longer reduce the +energy of the system. Individual atom forces may still be larger than +desired at this point, because the energy change is measured as the +difference of two large values (energy before and energy after) and +that difference may be smaller than machine epsilon even if atoms +could move in the gradient direction to reduce forces further.

By contast, the quadratic line search algorithm is often able to reduce forces closer to 0.0. It may also be more efficient than the diff -Naur lammps-28Sep09/doc/min_modify.txt lammps-29Sep09/doc/min_modify.txt --- lammps-28Sep09/doc/min_modify.txt 2009-02-25 17:18:18.000000000 -0700 +++ lammps-29Sep09/doc/min_modify.txt 2009-09-29 08:10:23.000000000 -0600 @@ -27,30 +27,29 @@ [Description:] This command sets parameters that affect the energy minimization -algorithms. The various settings may affect the convergence rate and -overall number of force evaluations required by a minimization, so -users can experiment with these parameters to tune their -minimizations. - -The minimization algorithms have an outer iteration (conjugate -gradient or steepest descent) and an inner iteration which is steps -along a one-dimensional line search in a particular search direction. -The {dmax} parameter is how far any atom can move in a single line -search in any dimension (x, y, or z). Thus a value of 0.1 in real -"units"_units.html means no atom will move further than 0.1 Angstroms -in a single outer iteration. This prevents highly overlapped atoms -from being moved long distances (e.g. through another atom) due to -large forces. - -The choice of line search algorithm can be selected via the {line} -keyword. The default backtracking search is very robust and should -always find a local energy minimum. However, it will "converge" when -it can no longer reduce the energy of the system. Individual atom -forces may still be larger than desired at this point, because the -energy change is measured as the difference of two large values -(energy before and energy after) and that difference may be smaller -than machine epsilon even if atoms could move in the gradient -direction to reduce forces further. +algorithms selected by the "min_style"_min_style.html command. The +various settings may affect the convergence rate and overall number of +force evaluations required by a minimization, so users can experiment +with these parameters to tune their minimizations. + +The {cg} and {sd} minimization styles have an outer iteration and an +inner iteration which is steps along a one-dimensional line search in +a particular search direction. The {dmax} parameter is how far any +atom can move in a single line search in any dimension (x, y, or z). +Thus a value of 0.1 in real "units"_units.html means no atom will move +further than 0.1 Angstroms in a single outer iteration. This prevents +highly overlapped atoms from being moved long distances (e.g. through +another atom) due to large forces. + +The choice of line search algorithm for the {cg} and {sd} minimization +styles can be selected via the {line} keyword. The default +backtracking search is robust and should always find a local energy +minimum. However, it will "converge" when it can no longer reduce the +energy of the system. Individual atom forces may still be larger than +desired at this point, because the energy change is measured as the +difference of two large values (energy before and energy after) and +that difference may be smaller than machine epsilon even if atoms +could move in the gradient direction to reduce forces further. By contast, the {quadratic} line search algorithm is often able to reduce forces closer to 0.0. It may also be more efficient than the diff -Naur lammps-28Sep09/doc/min_style.html lammps-29Sep09/doc/min_style.html --- lammps-28Sep09/doc/min_style.html 2008-04-14 15:37:10.000000000 -0600 +++ lammps-29Sep09/doc/min_style.html 2009-09-29 08:10:23.000000000 -0600 @@ -15,12 +15,12 @@ 

min_style style
-
  • style = cg or sd +
    • style = cg or hftn or sd

    Examples: 

    min_style cg
-min_style sd 
+min_style hftn

    Description:

    @@ -35,6 +35,18 @@ restarted when it ceases to make progress. The PR variant is thought to be the most effective CG choice.

    +

    Style hftn is a Hessian-free truncated Newton algorithm. At each +iteration a quadratic model of the energy potential is solved by a +conjugate gradient inner iteration. The Hessian (second derivatives) +of the energy is not formed directly, but approximated in each +conjugate search direction by a finite difference directional +derivative. When close to an energy minimum, the algorithm behaves +like a Newton method and exhibits a quadratic convergence rate to high +accuracy. In most cases the behavior of hftn is similar to cg, +but it offers another minimizer alternative if cg seems to perform +poorly. This style is not affected by the +min_modify command. +

    Style sd is a steepest descent algorithm. At each iteration, the search direction is set to the downhill direction corresponding to the force vector (negative gradient of energy). Typically, steepest diff -Naur lammps-28Sep09/doc/min_style.txt lammps-29Sep09/doc/min_style.txt --- lammps-28Sep09/doc/min_style.txt 2008-04-14 15:37:10.000000000 -0600 +++ lammps-29Sep09/doc/min_style.txt 2009-09-29 08:10:23.000000000 -0600 @@ -11,12 +11,12 @@ min_style style :pre -style = {cg} or {sd} :ul +style = {cg} or {hftn} or {sd} :ul [Examples:] min_style cg -min_style sd :pre +min_style hftn :pre [Description:] @@ -31,6 +31,18 @@ restarted when it ceases to make progress. The PR variant is thought to be the most effective CG choice. +Style {hftn} is a Hessian-free truncated Newton algorithm. At each +iteration a quadratic model of the energy potential is solved by a +conjugate gradient inner iteration. The Hessian (second derivatives) +of the energy is not formed directly, but approximated in each +conjugate search direction by a finite difference directional +derivative. When close to an energy minimum, the algorithm behaves +like a Newton method and exhibits a quadratic convergence rate to high +accuracy. In most cases the behavior of {hftn} is similar to {cg}, +but it offers another minimizer alternative if {cg} seems to perform +poorly. This style is not affected by the +"min_modify"_min_modify.html command. + Style {sd} is a steepest descent algorithm. At each iteration, the search direction is set to the downhill direction corresponding to the force vector (negative gradient of energy). Typically, steepest diff -Naur lammps-28Sep09/doc/shape.html lammps-29Sep09/doc/shape.html --- lammps-28Sep09/doc/shape.html 2009-06-26 12:22:33.000000000 -0600 +++ lammps-29Sep09/doc/shape.html 2009-09-29 08:10:23.000000000 -0600 @@ -85,8 +85,8 @@

    If you define a hybrid atom style which includes one (or more) sub-styles which require per-type shape and one (or more) sub-styles which require per-atom diameter, then you must define both. -In this case the per-type shape will be ignored; only the per-atom -diameter will be used by LAMMPS. Note that this means you can not +However, in this case the per-type shape will be ignored; only the +per-atom diameter will be used by LAMMPS. This means you cannot currently mix aspherical particles with per-atom diameter particles.

    Restrictions: diff -Naur lammps-28Sep09/doc/shape.txt lammps-29Sep09/doc/shape.txt --- lammps-28Sep09/doc/shape.txt 2009-06-26 12:22:33.000000000 -0600 +++ lammps-29Sep09/doc/shape.txt 2009-09-29 08:10:23.000000000 -0600 @@ -82,8 +82,8 @@ If you define a "hybrid atom style"_atom_style.html which includes one (or more) sub-styles which require per-type shape and one (or more) sub-styles which require per-atom diameter, then you must define both. -In this case the per-type shape will be ignored; only the per-atom -diameter will be used by LAMMPS. Note that this means you can not +However, in this case the per-type shape will be ignored; only the +per-atom diameter will be used by LAMMPS. This means you cannot currently mix aspherical particles with per-atom diameter particles. [Restrictions:] diff -Naur lammps-28Sep09/src/ASPHERE/fix_nve_asphere.cpp lammps-29Sep09/src/ASPHERE/fix_nve_asphere.cpp --- lammps-28Sep09/src/ASPHERE/fix_nve_asphere.cpp 2009-06-26 12:23:16.000000000 -0600 +++ lammps-29Sep09/src/ASPHERE/fix_nve_asphere.cpp 2009-09-29 08:10:33.000000000 -0600 @@ -27,9 +27,6 @@ #include "memory.h" #include "error.h" -#define TOLERANCE 1.0e-6 -#define EPSILON 1.0e-7 - using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ diff -Naur lammps-28Sep09/src/USER-ATC/README lammps-29Sep09/src/USER-ATC/README --- lammps-28Sep09/src/USER-ATC/README 1969-12-31 17:00:00.000000000 -0700 +++ lammps-29Sep09/src/USER-ATC/README 2009-09-28 15:28:19.000000000 -0600 @@ -0,0 +1,12 @@ +The files in this directory are a user-contributed package for LAMMPS. + +The primary people who created these files are Reese Jones +(rjones@sandia.gov), Jeremy Templeton (jatemple@sandia.gov) and Jon +Zimmerman (jzimmer@sandia.gov). Contact them directly if you have +questions. + +This package implements a "fix atc" command which can be used in a +LAMMPS input script. This fix can be employed to either do concurrent +coupling of MD with FE-based physics surrogates or on-the-fly +post-processing of atomic information to continuum fields. See the +doc page for the fix atc command to get started. diff -Naur lammps-28Sep09/src/fix.h lammps-29Sep09/src/fix.h --- lammps-28Sep09/src/fix.h 2009-08-28 14:28:23.000000000 -0600 +++ lammps-29Sep09/src/fix.h 2009-09-29 08:10:33.000000000 -0600 @@ -109,6 +109,9 @@ virtual double min_energy(double *) {return 0.0;} virtual void min_store() {} + virtual void min_clearstore() {} + virtual void min_pushstore() {} + virtual void min_popstore() {} virtual void min_step(double, double *) {} virtual double max_alpha(double *) {return 0.0;} virtual int min_dof() {return 0;} diff -Naur lammps-28Sep09/src/fix_box_relax.cpp lammps-29Sep09/src/fix_box_relax.cpp --- lammps-28Sep09/src/fix_box_relax.cpp 2009-08-18 11:56:26.000000000 -0600 +++ lammps-29Sep09/src/fix_box_relax.cpp 2009-09-29 08:10:33.000000000 -0600 @@ -179,6 +179,8 @@ dimension = domain->dimension; nrigid = 0; rfix = 0; + + current_lifo = 0; } /* ---------------------------------------------------------------------- */ @@ -301,20 +303,60 @@ } /* ---------------------------------------------------------------------- - store extra dof values for linesearch starting point + store extra dof values for minimization linesearch starting point boxlo0,boxhi0 = box dimensions s0 = ratio of current boxsize to initial boxsize + box values are pushed onto a LIFO stack so nested calls can be made + values are popped by calling min_step(0.0) ------------------------------------------------------------------------- */ void FixBoxRelax::min_store() { for (int i = 0; i < 3; i++) { - boxlo0[i] = domain->boxlo[i]; - boxhi0[i] = domain->boxhi[i]; + boxlo0[current_lifo][i] = domain->boxlo[i]; + boxhi0[current_lifo][i] = domain->boxhi[i]; + } + s0[0] = (boxhi0[current_lifo][0]-boxlo0[current_lifo][0])/xprdinit; + s0[1] = (boxhi0[current_lifo][1]-boxlo0[current_lifo][1])/yprdinit; + s0[2] = (boxhi0[current_lifo][2]-boxlo0[current_lifo][2])/zprdinit; +} + +/* ---------------------------------------------------------------------- + clear the LIFO stack for min_store +------------------------------------------------------------------------- */ + +void FixBoxRelax::min_clearstore() +{ + current_lifo = 0; +} + +/* ---------------------------------------------------------------------- + push the LIFO stack for min_store +------------------------------------------------------------------------- */ + +void FixBoxRelax::min_pushstore() +{ + if (current_lifo >= MAX_LIFO_DEPTH) { + error->all("Attempt to push beyond stack limit "); + return; + } + + current_lifo++; +} + + +/* ---------------------------------------------------------------------- + pop the LIFO stack for min_store +------------------------------------------------------------------------- */ + +void FixBoxRelax::min_popstore() +{ + if (current_lifo <= 0) { + error->all("Attempt to pop empty stack "); + return; } - s0[0] = (boxhi0[0]-boxlo0[0])/xprdinit; - s0[1] = (boxhi0[1]-boxlo0[1])/yprdinit; - s0[2] = (boxhi0[2]-boxlo0[2])/zprdinit; + + current_lifo--; } /* ---------------------------------------------------------------------- @@ -394,9 +436,11 @@ for (i = 0; i < 3; i++) if (p_flag[i]) { - ctr = 0.5 * (boxlo0[i] + boxhi0[i]); - domain->boxlo[i] = boxlo0[i] + (boxlo0[i]-ctr)*ds[i]/s0[i]; - domain->boxhi[i] = boxhi0[i] + (boxhi0[i]-ctr)*ds[i]/s0[i]; + double currentBoxLo0 = boxlo0[current_lifo][i]; + double currentBoxHi0 = boxhi0[current_lifo][i]; + ctr = 0.5 * (currentBoxLo0 + currentBoxHi0); + domain->boxlo[i] = currentBoxLo0 + (currentBoxLo0-ctr)*ds[i]/s0[i]; + domain->boxhi[i] = currentBoxHi0 + (currentBoxHi0-ctr)*ds[i]/s0[i]; } domain->set_global_box(); diff -Naur lammps-28Sep09/src/fix_box_relax.h lammps-29Sep09/src/fix_box_relax.h --- lammps-28Sep09/src/fix_box_relax.h 2009-08-14 14:54:10.000000000 -0600 +++ lammps-29Sep09/src/fix_box_relax.h 2009-09-29 08:10:33.000000000 -0600 @@ -16,6 +16,8 @@ #include "fix.h" +const int MAX_LIFO_DEPTH = 2; // to meet the needs of min_hftn + namespace LAMMPS_NS { class FixBoxRelax : public Fix { @@ -27,6 +29,9 @@ double min_energy(double *); void min_store(); + void min_clearstore(); + void min_pushstore(); + void min_popstore(); void min_step(double, double *); double max_alpha(double *); int min_dof(); @@ -42,9 +47,11 @@ double volinit,xprdinit,yprdinit,zprdinit; double vmax,pv2e,pflagsum; - double boxlo0[3],boxhi0[3]; // box bounds at start of line search - double s0[6]; // scale matrix in Voigt notation - double ds[6]; // increment in scale matrix + int current_lifo; // LIFO stack pointer + double boxlo0[MAX_LIFO_DEPTH][3]; // box bounds at start of line search + double boxhi0[MAX_LIFO_DEPTH][3]; + double s0[6]; // scale matrix in Voigt notation + double ds[6]; // increment in scale matrix char *id_temp,*id_press; class Compute *temperature,*pressure; diff -Naur lammps-28Sep09/src/min.cpp lammps-29Sep09/src/min.cpp --- lammps-28Sep09/src/min.cpp 2009-08-28 14:28:23.000000000 -0600 +++ lammps-29Sep09/src/min.cpp 2009-09-29 08:10:33.000000000 -0600 @@ -304,11 +304,16 @@ { // possible stop conditions - char *stopstrings[] = {"max iterations","max force evaluations", - "energy tolerance","force tolerance", - "search direction is not downhill", - "linesearch alpha is zero", - "forces are zero","quadratic factors are zero"}; + char *stopstrings[] = {"max iterations", + "max force evaluations", + "energy tolerance", + "force tolerance", + "search direction is not downhill", + "linesearch alpha is zero", + "forces are zero", + "quadratic factors are zero", + "trust region too small", + "HFTN minimizer error"}; // stats for Finish to print diff -Naur lammps-28Sep09/src/min_hftn.cpp lammps-29Sep09/src/min_hftn.cpp --- lammps-28Sep09/src/min_hftn.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-29Sep09/src/min_hftn.cpp 2009-09-29 08:10:33.000000000 -0600 @@ -0,0 +1,1680 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Author: Todd Plantenga (SNL) + Sources: "Numerical Optimization", Nocedal and Wright, 2nd Ed, p170 + "Parallel Unconstrained Min", Plantenga, SAND98-8201 +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "atom.h" +#include "fix_minimize.h" +#include "min_hftn.h" +#include "modify.h" +#include "output.h" +#include "pair.h" +#include "update.h" +#include "timer.h" + +#define MIN(A,B) (((A) < (B)) ? (A) : (B)) +#define MAX(A,B) (((A) > (B)) ? (A) : (B)) + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + * This class performs Hessian-free truncated Newton minimization on an + * unconstrained molecular potential. The algorithm avoids computing the + * Hessian matrix, but obtains a near-quadratic rate of convergence. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + File local data +------------------------------------------------------------------------- */ + +//---- CONSTANTS MAP TO stopstrings DECLARED IN Min.run (min.cpp). + +static const int STOP_MAX_ITER = 0; //-- MAX ITERATIONS EXCEEDED +static const int STOP_MAX_FORCE_EVALS = 1; //-- MAX FORCE EVALUATIONS EXCEEDED +static const int STOP_ENERGY_TOL = 2; //-- STEP DID NOT CHANGE ENERGY +static const int STOP_FORCE_TOL = 3; //-- CONVERGED TO DESIRED FORCE TOL +static const int STOP_TR_TOO_SMALL = 8; //-- TRUST REGION TOO SMALL +static const int STOP_ERROR = 9; //-- INTERNAL ERROR + +static const int NO_CGSTEP_BECAUSE_F_TOL_SATISFIED = 0; +static const int CGSTEP_NEWTON = 1; +static const int CGSTEP_TO_TR = 2; +static const int CGSTEP_TO_DMAX = 3; +static const int CGSTEP_NEGATIVE_CURVATURE = 4; +static const int CGSTEP_MAX_INNER_ITERS = 5; +static const int CGSTEP_UNDETERMINED = 6; + +//---- WHEN TESTING ENERGY_TOL, THE ENERGY MAGNITUDE MUST BE AT LEAST THIS BIG. + +static const double MIN_ETOL_MAG = 1.0e-8; + +//---- MACHINE PRECISION IS SOMETIMES DEFINED BY THE C RUNTIME. + +#ifdef DBL_EPSILON + #define MACHINE_EPS DBL_EPSILON +#else + #define MACHINE_EPS 2.220446049250313e-16 +#endif + +/* ---------------------------------------------------------------------- + Constructor +------------------------------------------------------------------------- */ + +MinHFTN::MinHFTN(LAMMPS *lmp) : Min(lmp) +{ + for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) + _daExtraGlobal[i] = NULL; + for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) + _daExtraAtom[i] = NULL; + + _fpPrint = NULL; + + return; +} + +/* ---------------------------------------------------------------------- + Destructor +------------------------------------------------------------------------- */ + +MinHFTN::~MinHFTN (void) +{ + for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) + if (_daExtraGlobal[i] != NULL) + delete [] _daExtraGlobal[i]; + for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) + if (_daExtraAtom[i] != NULL) + delete [] _daExtraAtom[i]; + + return; +} + +/* ---------------------------------------------------------------------- + Public method init_style +------------------------------------------------------------------------- */ + +void MinHFTN::init_style() +{ + for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) { + if (_daExtraGlobal[i] != NULL) + delete [] _daExtraGlobal[i]; + _daExtraGlobal[i] = NULL; + } + for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) { + if (_daExtraAtom[i] != NULL) + delete [] _daExtraAtom[i]; + _daExtraAtom[i] = NULL; + } + + return; +} + +/* ---------------------------------------------------------------------- + Public method setup_style +------------------------------------------------------------------------- */ + +void MinHFTN::setup_style() +{ + //---- ALLOCATE MEMORY FOR ATOMIC DEGREES OF FREEDOM. + for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) + fix_minimize->add_vector(3); + + //---- ALLOCATE MEMORY FOR EXTRA GLOBAL DEGREES OF FREEDOM. + //---- THE FIX MODULE TAKES CARE OF THE FIRST VECTOR, X0 (XK). + if (nextra_global) { + for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) + _daExtraGlobal[i] = new double[nextra_global]; + } + + //---- ALLOCATE MEMORY FOR EXTRA PER-ATOM DEGREES OF FREEDOM. + if (nextra_atom) { + for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) + _daExtraAtom[i] = new double*[nextra_atom]; + + for (int m = 0; m < nextra_atom; m++) { + for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) + fix_minimize->add_vector (extra_peratom[m]); + } + } + + return; +} + +/* ---------------------------------------------------------------------- + Public method reset_vectors + After an energy/force calculation, atoms may migrate from one processor + to another. Any local vector correlated with atom positions or forces + must also be migrated. This is accomplished by a subclass of Fix. + This method updates local pointers to the latest Fix copies. +------------------------------------------------------------------------- */ + +void MinHFTN::reset_vectors() +{ + n3 = 3 * atom->nlocal; + + //---- ATOMIC DEGREES OF FREEDOM. + if (n3 > 0) { + x = atom->x[0]; + f = atom->f[0]; + } + for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) + _daAVectors[i] = fix_minimize->request_vector (i); + + //---- EXTRA PER-ATOM DEGREES OF FREEDOM. + if (nextra_atom) { + int n = NUM_HFTN_ATOM_BASED_VECTORS; + for (int m = 0; m < nextra_atom; m++) { + extra_nlen[m] = extra_peratom[m] * atom->nlocal; + requestor[m]->min_pointers (&xextra_atom[m], &fextra_atom[m]); + for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) + _daExtraAtom[i][m] = fix_minimize->request_vector (n++); + } + } + + return; +} + +/* ---------------------------------------------------------------------- + Public method iterate + Upon entry, Min::setup() and Min::run have executed, and energy has + already been evaluated at the initial point. Return an integer code + that maps to a stop condition in min.cpp. +------------------------------------------------------------------------- */ + +int MinHFTN::iterate(int) +{ + //---- TURN THIS ON TO GENERATE AN OPTIMIZATION PROGRESS FILE. + bool bPrintProgress = false; + + if (bPrintProgress) + open_hftn_print_file_(); + + double dFinalEnergy = 0.0; + double dFinalFnorm2 = 0.0; + modify->min_clearstore(); + int nStopCode = execute_hftn_ (bPrintProgress, + einitial, + fnorm2_init, + dFinalEnergy, + dFinalFnorm2); + modify->min_clearstore(); + if (bPrintProgress) + close_hftn_print_file_(); + + return( nStopCode ); +} + +/* ---------------------------------------------------------------------- + Private method execute_hftn_ + @param[in] bPrintProgress - if true then print progress to a file + @param[in] dInitialEnergy - energy at input x + @param[in] dInitialForce2 - |F|_2 at input x + @param[out] dFinalEnergy - energy at output x + @param[out] dFinalForce2 - |F|_2 at output x + + Return stop code described in the enumeration at the top of this file, + and the following: + atom->x - positions at output x + atom->f - forces evaluated at output x +------------------------------------------------------------------------- */ +int MinHFTN::execute_hftn_(const bool bPrintProgress, + const double dInitialEnergy, + const double dInitialForce2, + double & dFinalEnergy, + double & dFinalForce2) +{ + //---- DEFINE OUTPUTS PRINTED BY "Finish". + eprevious = dInitialEnergy; + alpha_final = 0.0; + neval = 0; + niter = 0; + dFinalEnergy = dInitialEnergy; + dFinalForce2 = dInitialForce2; + + if (dInitialForce2 < update->ftol) + return( STOP_FORCE_TOL ); + + //---- SAVE ATOM POSITIONS BEFORE AN ITERATION. + fix_minimize->store_box(); + for (int i = 0; i < n3; i++) + _daAVectors[VEC_XK][i] = atom->x[0][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * xkAtom = _daExtraAtom[VEC_XK][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xkAtom[i] = xatom[i]; + } + } + if (nextra_global) + modify->min_store(); + + double dXInf = calc_xinf_using_mpi_(); + + //---- FIND THE NUMBER OF UNKNOWNS. + int nLocalNumUnknowns = n3 + nextra_atom; + MPI_Allreduce (&nLocalNumUnknowns, &_nNumUnknowns, + 1, MPI_INT, MPI_SUM, world); + + //---- INITIALIZE THE TRUST RADIUS BASED ON THE GRADIENT. + double dTrustRadius = 1.5 * dInitialForce2; + + //---- TRUST RADIUS MUST KEEP STEPS FROM LETTING ATOMS MOVE SO FAR THEY + //---- VIOLATE PHYSICS OR JUMP BEYOND A PARALLEL PROCESSING DOMAIN. + //---- LINE SEARCH METHODS DO THIS BY RESTRICTING THE LARGEST CHANGE + //---- OF ANY ATOM'S COMPONENT TO dmax. AN EXACT CHECK IS MADE LATER, + //---- BUT THIS GUIDES DETERMINATION OF A MAX TRUST RADIUS. + double dMaxTrustRadius = dmax * sqrt (_nNumUnknowns); + + dTrustRadius = MIN (dTrustRadius, dMaxTrustRadius); + double dLastNewtonStep2 = dMaxTrustRadius; + + if (bPrintProgress) + hftn_print_line_ (false, -1, neval, dInitialEnergy, dInitialForce2, + -1, dTrustRadius, 0.0, 0.0, 0.0); + + bool bHaveEvaluatedAtX = true; + double dCurrentEnergy = dInitialEnergy; + double dCurrentForce2 = dInitialForce2; + for (niter = 0; niter < update->nsteps; niter++) { + (update->ntimestep)++; + + //---- CALL THE INNER LOOP TO GET THE NEXT TRUST REGION STEP. + + double dCgForce2StopTol = MIN ((dCurrentForce2 / 2.0), 0.1 / (niter+1)); + dCgForce2StopTol = MAX (dCgForce2StopTol, update->ftol); + + double dNewEnergy; + double dNewForce2; + int nStepType; + double dStepLength2; + double dStepLengthInf; + if (compute_inner_cg_step_ (dTrustRadius, + dCgForce2StopTol, + update->max_eval, + bHaveEvaluatedAtX, + dCurrentEnergy, dCurrentForce2, + dNewEnergy, dNewForce2, + nStepType, + dStepLength2, dStepLengthInf) == false) { + //---- THERE WAS AN ERROR. RESTORE TO LAST ACCEPTED STEP. + if (nextra_global) + modify->min_step (0.0, _daExtraGlobal[VEC_CG_P]); + for (int i = 0; i < n3; i++) + atom->x[0][i] = _daAVectors[VEC_XK][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * xkAtom = _daExtraAtom[VEC_XK][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xatom[i] = xkAtom[i]; + } + } + dFinalEnergy = energy_force (0); + neval++; + dFinalForce2 = sqrt (fnorm_sqr()); + return( STOP_ERROR ); + } + + //---- STOP IF THE CURRENT POSITION WAS FOUND TO BE ALREADY GOOD ENOUGH. + //---- IN THIS CASE THE ENERGY AND FORCES ARE ALREADY COMPUTED. + if (nStepType == NO_CGSTEP_BECAUSE_F_TOL_SATISFIED) { + if (bPrintProgress) + hftn_print_line_ (true, niter+1, neval, dNewEnergy, dNewForce2, + nStepType, dTrustRadius, dStepLength2, + 0.0, 0.0); + dFinalEnergy = dNewEnergy; + dFinalForce2 = dNewForce2; + return( STOP_FORCE_TOL ); + } + + //---- COMPUTE THE DIRECTIONAL DERIVATIVE H(x_k) p. + bool bUseForwardDiffs = (dCurrentForce2 > 1000.0 * sqrt (MACHINE_EPS)); + evaluate_dir_der_ (bUseForwardDiffs, + VEC_CG_P, + VEC_CG_HD, + true, + dCurrentEnergy); + + //---- COMPUTE p^T grad(x_k) AND SAVE IT FOR PRED. + double dGradDotP = calc_grad_dot_v_using_mpi_ (VEC_CG_P); + + //---- MOVE TO THE NEW POINT AND EVALUATE ENERGY AND FORCES. + //---- THIS IS THE PLACE WHERE energy_force IS ALLOWED TO RESET. + for (int i = 0; i < n3; i++) + atom->x[0][i] = _daAVectors[VEC_XK][i] + _daAVectors[VEC_CG_P][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * xkAtom = _daExtraAtom[VEC_XK][m]; + double * pAtom = _daExtraAtom[VEC_CG_P][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xatom[i] = xkAtom[i] + pAtom[i]; + } + } + if (nextra_global) + modify->min_step (1.0, _daExtraGlobal[VEC_CG_P]); + dNewEnergy = energy_force (1); + neval++; + + dNewForce2 = sqrt (fnorm_sqr()); + + double dAred = dCurrentEnergy - dNewEnergy; + + //---- STOP IF THE FORCE TOLERANCE IS MET. + if (dNewForce2 < update->ftol) { + if (bPrintProgress) + hftn_print_line_ (true, niter+1, neval, dNewEnergy, dNewForce2, + nStepType, dTrustRadius, dStepLength2, + dAred, -1.0); + //---- (IMPLICITLY ACCEPT THE LAST STEP TO THE NEW POINT.) + dFinalEnergy = dNewEnergy; + dFinalForce2 = dNewForce2; + return( STOP_FORCE_TOL ); + } + + //---- STOP IF THE ACTUAL ENERGY REDUCTION IS TINY. + if (nStepType != CGSTEP_TO_DMAX) { + double dMag = 0.5 * (fabs (dCurrentEnergy) + fabs (dNewEnergy)); + dMag = MAX (dMag, MIN_ETOL_MAG); + if ( (fabs (dAred) < (update->etol * dMag)) + || (dStepLengthInf == 0.0) ) { + if (bPrintProgress) + hftn_print_line_ (true, niter+1, neval, + dNewEnergy, dNewForce2, + nStepType, dTrustRadius, dStepLength2, + dAred, -1.0); + //---- (IMPLICITLY ACCEPT THE LAST STEP TO THE NEW POINT.) + dFinalEnergy = dNewEnergy; + dFinalForce2 = dNewForce2; + return( STOP_ENERGY_TOL ); + } + } + + //---- COMPUTE THE PREDICTED REDUCTION - p^T grad - 0.5 p^T Hp + double dPHP = calc_dot_prod_using_mpi_ (VEC_CG_P, VEC_CG_HD); + double dPred = - dGradDotP - (0.5 * dPHP); + + //---- ACCEPT OR REJECT THE STEP PROPOSED BY THE INNER CG LOOP. + //---- WHEN NEAR A SOLUTION, THE FORCE NORM IS PROBABLY MORE ACCURATE, + //---- SO DON'T ACCEPT A STEP THAT REDUCES ENERGY SOME TINY AMOUNT + //---- WHILE INCREASING THE FORCE NORM. + bool bStepAccepted = (dAred > 0.0) + && ( (dNewForce2 < dCurrentForce2) + || (dCurrentForce2 > 1.0e-6)); + if (bStepAccepted) { + //---- THE STEP IS ACCEPTED. + if (bPrintProgress) + hftn_print_line_ (true, niter+1, neval, dNewEnergy, dNewForce2, + nStepType, dTrustRadius, dStepLength2, + dAred, dPred); + + fix_minimize->store_box(); + modify->min_clearstore(); + for (int i = 0; i < n3; i++) + _daAVectors[VEC_XK][i] = atom->x[0][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * xkAtom = _daExtraAtom[VEC_XK][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xkAtom[i] = xatom[i]; + } + } + if (nextra_global) + modify->min_store(); + + if (niter > 0) + eprevious = dCurrentEnergy; + dCurrentEnergy = dNewEnergy; + dCurrentForce2 = dNewForce2; + bHaveEvaluatedAtX = true; + + if (nStepType == CGSTEP_NEWTON) + dLastNewtonStep2 = dStepLength2; + + //---- UPDATE THE TRUST REGION BASED ON AGREEMENT BETWEEN + //---- THE ACTUAL REDUCTION AND THE PREDICTED MODEL REDUCTION. + if ((dAred > 0.75 * dPred) && (dStepLength2 >= 0.99 * dTrustRadius)) + dTrustRadius = 2.0 * dTrustRadius; + dTrustRadius = MIN (dTrustRadius, dMaxTrustRadius); + + //---- DMAX VIOLATIONS TRUNCATE THE CG STEP WITHOUT COMPARISONS; + //---- BETTER TO ADJUST THE TRUST REGION SO DMAX STOPS HAPPENING. + if (nStepType == CGSTEP_TO_DMAX) { + if (dStepLength2 <= MACHINE_EPS) + dTrustRadius = 0.1 * dTrustRadius; + else + dTrustRadius = MIN (dTrustRadius, 2.0 * dStepLength2); + } + } + else { + //---- THE STEP IS REJECTED. + if (bPrintProgress) + hftn_print_line_ (false, niter+1, neval, + dCurrentEnergy, dCurrentForce2, + nStepType, dTrustRadius, dStepLength2, + dAred, dPred); + + //---- RESTORE THE LAST X_K POSITION. + if (nextra_global) + modify->min_step (0.0, _daExtraGlobal[VEC_CG_P]); + for (int i = 0; i < n3; i++) + atom->x[0][i] = _daAVectors[VEC_XK][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * xkAtom = _daExtraAtom[VEC_XK][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xatom[i] = xkAtom[i]; + } + } + bHaveEvaluatedAtX = false; + + //---- UPDATE THE TRUST REGION. + //---- EXPERIMENTS INDICATE NEGATIVE CURVATURE CAN TAKE A BAD + //---- STEP A LONG WAY, SO BE MORE AGGRESSIVE IN THIS CASE. + //---- ALSO, IF NEAR A SOLUTION AND DONE WITH NEWTON STEPS, + //---- THEN REDUCE TO SOMETHING NEAR THE LAST GOOD NEWTON STEP. + if ((nStepType == CGSTEP_NEGATIVE_CURVATURE) && (-dAred > dPred)) + dTrustRadius = 0.10 * MIN (dTrustRadius, dStepLength2); + else if ( (nStepType == CGSTEP_TO_DMAX) + && (dStepLength2 <= MACHINE_EPS)) + dTrustRadius = 0.10 * dTrustRadius; + else if (-dAred > dPred) + dTrustRadius = 0.20 * MIN (dTrustRadius, dStepLength2); + else + dTrustRadius = 0.25 * MIN (dTrustRadius, dStepLength2); + + if ( (nStepType != CGSTEP_NEWTON) + && (dCurrentForce2 < sqrt (MACHINE_EPS))) + dTrustRadius = MIN (dTrustRadius, 2.0 * dLastNewtonStep2); + + dLastNewtonStep2 = dMaxTrustRadius; + + //---- STOP IF THE TRUST RADIUS IS TOO SMALL TO CONTINUE. + if ( (dTrustRadius <= 0.0) + || (dTrustRadius <= MACHINE_EPS * MAX (1.0, dXInf))) { + dFinalEnergy = dCurrentEnergy; + dFinalForce2 = dCurrentForce2; + return( STOP_TR_TOO_SMALL ); + } + } + + //---- OUTPUT FOR thermo, dump, restart FILES. + if (output->next == update->ntimestep) { + //---- IF THE LAST STEP WAS REJECTED, THEN REEVALUATE ENERGY AND + //---- FORCES AT THE OLD POINT SO THE OUTPUT DOES NOT DISPLAY + //---- THE INCREASED ENERGY OF THE REJECTED STEP. + if (bStepAccepted == false) { + dCurrentEnergy = energy_force (1); + neval++; + } + timer->stamp(); + output->write (update->ntimestep); + timer->stamp (TIME_OUTPUT); + } + + //---- RETURN IF NUMBER OF EVALUATIONS EXCEEDED. + if (neval >= update->max_eval) { + dFinalEnergy = dCurrentEnergy; + dFinalForce2 = dCurrentForce2; + return( STOP_MAX_FORCE_EVALS ); + } + + } //-- END for LOOP OVER niter + + dFinalEnergy = dCurrentEnergy; + dFinalForce2 = dCurrentForce2; + return( STOP_MAX_ITER ); +} + +/* ---------------------------------------------------------------------- + Private method compute_inner_cg_step_ + Execute CG using Hessian-vector products approximated by finite difference + directional derivatives. + + On input these must be defined: + atom->x - positions at x + atom->f - ignored + VEC_XK - positions at x + On output these are defined: + atom->x - unchanged + atom->f - forces evaluated at x, but only if nStepType == NO_CGSTEP + VEC_XK - unchanged + VEC_CG_P - step from VEC_XK to new positions + During processing these are modified: + VEC_CG_D - conjugate gradient inner loop step + VEC_CG_HD - Hessian-vector product + VEC_CG_R - residual of inner loop step + VEC_DIF1 - temp storage + VEC_DIF2 - temp storage + + @param[in] dTrustRadius - trust region radius for this subiteration + @param[in] dForceTol - stop tolerance on |F|_2 for this subiteration + @param[in] nMaxEvals - total energy/force evaluations allowed + @param[in] bHaveEvalAtXin - true if forces are valid at input x + @param[in] dEnergyAtXin - energy at input x, if bHaveEvalAtXin is true + @param[in] dForce2AtXin - |F|_2 at input x, if bHaveEvalAtXin is true + @param[out] dEnergyAtXout - energy at output x, if NO_CGSTEP (see below) + @param[out] dForce2AtXout - |F|_2 at output x, if NO_CGSTEP (see below) + @param[out] nStepType - step type for hftn_print_line_() + @param[out] dStepLength2 - |step|_2 + @param[out] dStepLengthInf - |step|_inf + + Return false if there was a fatal error. + If nStepType equals NO_CGSTEP_BECAUSE_F_TOL_SATISFIED, then the energy + and forces are evaluated and returned in dEnergyAtXout, dForce2AtXout; + else energy and forces are not evaluated. +------------------------------------------------------------------------- */ + +bool MinHFTN::compute_inner_cg_step_(const double dTrustRadius, + const double dForceTol, + const int nMaxEvals, + const bool bHaveEvalAtXin, + const double dEnergyAtXin, + const double dForce2AtXin, + double & dEnergyAtXout, + double & dForce2AtXout, + int & nStepType, + double & dStepLength2, + double & dStepLengthInf) +{ + //---- SET p_0 = 0. + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) + _daExtraGlobal[VEC_CG_P][i] = 0.0; + } + for (int i = 0; i < n3; i++) + _daAVectors[VEC_CG_P][i] = 0.0; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * pAtom = _daExtraAtom[VEC_CG_P][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + pAtom[i] = 0.0; + } + } + double dPP = 0.0; + + //---- OBTAIN THE ENERGY AND FORCES AT THE INPUT POSITION. + double dEnergyAtX = dEnergyAtXin; + double dForce2AtX = dForce2AtXin; + if (bHaveEvalAtXin == false) { + dEnergyAtX = energy_force (0); + neval++; + dForce2AtX = sqrt (fnorm_sqr()); + } + + //---- RETURN IMMEDIATELY IF THE FORCE TOLERANCE IS ALREADY MET. + //---- THE STEP TYPE INFORMS THE CALLER THAT ENERGY AND FORCES HAVE + //---- BEEN EVALUATED. + if (dForce2AtX <= dForceTol) { + dEnergyAtXout = dEnergyAtX; + dForce2AtXout = dForce2AtX; + nStepType = NO_CGSTEP_BECAUSE_F_TOL_SATISFIED; + dStepLength2 = 0.0; + dStepLengthInf = 0.0; + return( true ); + } + + //---- r_0 = -grad (FIRST SEARCH DIRECTION IS STEEPEST DESCENT) + //---- d_0 = r_0 + //---- REMEMBER THAT FORCES = -GRADIENT. + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) { + _daExtraGlobal[VEC_CG_R][i] = fextra[i]; + _daExtraGlobal[VEC_CG_D][i] = fextra[i]; + } + } + for (int i = 0; i < n3; i++) { + _daAVectors[VEC_CG_R][i] = atom->f[0][i]; + _daAVectors[VEC_CG_D][i] = atom->f[0][i]; + } + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * rAtom = _daExtraAtom[VEC_CG_R][m]; + double * dAtom = _daExtraAtom[VEC_CG_D][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) { + rAtom[i] = fextra_atom[m][i]; + dAtom[i] = fextra_atom[m][i]; + } + } + } + double dRR = dForce2AtX * dForce2AtX; + double dR0norm2 = sqrt (dRR); + + //---- LIMIT THE NUMBER OF INNER CG ITERATIONS. + //---- BASE IT ON THE NUMBER OF UNKNOWNS, OR MAXIMUM EVALUATIONS ASSUMING + //---- FORWARD DIFFERENCES ARE USED. + //---- NOTE THAT SETTING MAX=1 GIVES STEEPEST DESCENT. + int nLimit1 = _nNumUnknowns / 5; + if (nLimit1 < 100) + nLimit1 = MIN (_nNumUnknowns, 100); + int nLimit2 = (nMaxEvals - neval) / 2; + int nMaxInnerIters = MIN (nLimit1, nLimit2); + + //---- FURTHER LIMIT ITERATIONS IF NEAR MACHINE ROUNDOFF. + //---- THE METHOD CAN WASTE A LOT EVALUATIONS WITH LITTLE PAYOFF PROSPECT. + if (dForce2AtX < (sqrt (MACHINE_EPS) * MAX (1.0, fabs (dEnergyAtX))) ) + nMaxInnerIters = MIN (nMaxInnerIters, _nNumUnknowns / 20); + + bool bUseForwardDiffs = (dForce2AtX > 1000.0 * sqrt (MACHINE_EPS)); + + //---- MAIN CG LOOP. + for (int nInnerIter = 0; nInnerIter < nMaxInnerIters; nInnerIter++) { + //---- COMPUTE HESSIAN-VECTOR PRODUCT: H(x_k) d_i. + double dDummyEnergy; + evaluate_dir_der_ (bUseForwardDiffs, + VEC_CG_D, + VEC_CG_HD, + false, + dDummyEnergy); + + //---- CALCULATE d_i^T H d_i AND d_i^T d_i. + double dDHD; + double dDD; + calc_dhd_dd_using_mpi_ (dDHD, dDD); + + //---- HANDLE NEGATIVE CURVATURE. + if (dDHD <= (MACHINE_EPS * dDD)) { + //---- PROJECT BOTH DIRECTIONS TO THE TRUST RADIUS AND DECIDE + //---- WHICH MAKES A BETTER PREDICTED REDUCTION. + //---- p_i^T H(x_k) d_i AND grad_i^T d_i. + + double dPdotD = calc_dot_prod_using_mpi_ (VEC_CG_P, VEC_CG_D); + double dPdotHD = calc_dot_prod_using_mpi_ (VEC_CG_P, VEC_CG_HD); + + //---- MOVE TO X_K AND COMPUTE ENERGY AND FORCES. + if (nextra_global) + modify->min_step (0.0, _daExtraGlobal[VEC_CG_P]); + for (int i = 0; i < n3; i++) + atom->x[0][i] = _daAVectors[VEC_XK][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * xkAtom = _daExtraAtom[VEC_XK][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xatom[i] = xkAtom[i]; + } + } + dEnergyAtX = energy_force (0); + neval++; + + double dGradDotD = calc_grad_dot_v_using_mpi_ (VEC_CG_D); + + double tau = compute_to_tr_ (dPP, dPdotD, dDD, dTrustRadius, + true, dDHD, dPdotHD, dGradDotD); + + //---- MOVE THE POINT. + if (nextra_global) { + double * pGlobal = _daExtraGlobal[VEC_CG_P]; + double * dGlobal = _daExtraGlobal[VEC_CG_D]; + for (int i = 0; i < nextra_global; i++) { + pGlobal[i] += tau * dGlobal[i]; + } + } + for (int i = 0; i < n3; i++) + _daAVectors[VEC_CG_P][i] += tau * _daAVectors[VEC_CG_D][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * pAtom = _daExtraAtom[VEC_CG_P][m]; + double * dAtom = _daExtraAtom[VEC_CG_D][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + pAtom[i] += tau * dAtom[i]; + } + } + + nStepType = CGSTEP_NEGATIVE_CURVATURE; + calc_plengths_using_mpi_ (dStepLength2, dStepLengthInf); + return( true ); + } + + //---- COMPUTE THE OPTIMAL STEP LENGTH BASED ON THE QUADRATIC CG MODEL. + double dAlpha = dRR / dDHD; + + //---- MIGHT WANT TO ENABLE THIS TO DEBUG INTERNAL CG STEPS. + //fprintf (_fpPrint, " alpha = %11.8f neval=%4d\n", dAlpha, neval); + + //---- p_i+1 = p_i + alpha_i d_i + //---- (SAVE THE CURRENT p_i IN CASE THE STEP HAS TO BE SHORTENED.) + if (nextra_global) { + double * pGlobal = _daExtraGlobal[VEC_CG_P]; + double * dGlobal = _daExtraGlobal[VEC_CG_D]; + double * d1Global = _daExtraGlobal[VEC_DIF1]; + for (int i = 0; i < nextra_global; i++) { + d1Global[i] = pGlobal[i]; + pGlobal[i] += dAlpha * dGlobal[i]; + } + } + for (int i = 0; i < n3; i++) { + _daAVectors[VEC_DIF1][i] = _daAVectors[VEC_CG_P][i]; + _daAVectors[VEC_CG_P][i] += dAlpha * _daAVectors[VEC_CG_D][i]; + } + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * pAtom = _daExtraAtom[VEC_CG_P][m]; + double * dAtom = _daExtraAtom[VEC_CG_D][m]; + double * d1Atom = _daExtraAtom[VEC_DIF1][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) { + d1Atom[i] = pAtom[i]; + pAtom[i] += dAlpha * dAtom[i]; + } + } + } + + //---- COMPUTE VECTOR PRODUCTS p_i+1^T p_i+1 AND p_i^T d_i. + double dPnewDotPnew; + double dPoldDotD; + calc_ppnew_pdold_using_mpi_ (dPnewDotPnew, dPoldDotD); + + nStepType = CGSTEP_UNDETERMINED; + + //---- IF STEP LENGTH IS TOO LARGE, THEN REDUCE IT AND RETURN. + double tau; + if (step_exceeds_TR_ (dTrustRadius, dPP, dPoldDotD, dDD, tau)) { + adjust_step_to_tau_ (tau); + nStepType = CGSTEP_TO_TR; + } + if (step_exceeds_DMAX_()) { + adjust_step_to_tau_ (0.0); + nStepType = CGSTEP_TO_DMAX; + } + if ((nStepType == CGSTEP_TO_TR) || (nStepType == CGSTEP_TO_DMAX)) { + calc_plengths_using_mpi_ (dStepLength2, dStepLengthInf); + return( true ); + } + + dStepLength2 = sqrt (dPnewDotPnew); + + //---- r_i+1 = r_i - alpha * H d_i + if (nextra_global) { + double * rGlobal = _daExtraGlobal[VEC_CG_R]; + double * hdGlobal = _daExtraGlobal[VEC_CG_HD]; + for (int i = 0; i < nextra_global; i++) + rGlobal[i] -= dAlpha * hdGlobal[i]; + } + for (int i = 0; i < n3; i++) + _daAVectors[VEC_CG_R][i] -= dAlpha * _daAVectors[VEC_CG_HD][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * rAtom = _daExtraAtom[VEC_CG_R][m]; + double * hdAtom = _daExtraAtom[VEC_CG_HD][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + rAtom[i] -= dAlpha * hdAtom[i]; + } + } + double dRnewDotRnew = calc_dot_prod_using_mpi_ (VEC_CG_R, VEC_CG_R); + + //---- IF RESIDUAL IS SMALL ENOUGH, THEN RETURN THE CURRENT STEP. + if (sqrt (dRnewDotRnew) < dForceTol * dR0norm2) { + nStepType = CGSTEP_NEWTON; + calc_plengths_using_mpi_ (dStepLength2, dStepLengthInf); + return( true ); + } + + //---- beta = r_i+1^T r_i+1 / r_i^T r_i + //---- d_i+1 = r_i+1 + beta d_i + double dBeta = dRnewDotRnew / dRR; + if (nextra_global) { + double * rGlobal = _daExtraGlobal[VEC_CG_R]; + double * dGlobal = _daExtraGlobal[VEC_CG_D]; + for (int i = 0; i < nextra_global; i++) + dGlobal[i] = rGlobal[i] + dBeta * dGlobal[i]; + } + for (int i = 0; i < n3; i++) + _daAVectors[VEC_CG_D][i] = _daAVectors[VEC_CG_R][i] + + dBeta * _daAVectors[VEC_CG_D][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * rAtom = _daExtraAtom[VEC_CG_R][m]; + double * dAtom = _daExtraAtom[VEC_CG_D][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + dAtom[i] = rAtom[i] + dBeta * dAtom[i]; + } + } + + //---- CONTINUE THE LOOP. + dRR = dRnewDotRnew; + dPP = dPnewDotPnew; + } + + nStepType = CGSTEP_MAX_INNER_ITERS; + calc_plengths_using_mpi_ (dStepLength2, dStepLengthInf); + return( true ); +} + +/* ---------------------------------------------------------------------- + Private method calc_xinf_using_mpi_ +------------------------------------------------------------------------- */ + +double MinHFTN::calc_xinf_using_mpi_(void) const +{ + double dXInfLocal = 0.0; + for (int i = 0; i < n3; i++) + dXInfLocal = MAX (dXInfLocal, fabs (atom->x[0][i])); + + double dXInf; + MPI_Allreduce (&dXInfLocal, &dXInf, 1, MPI_DOUBLE, MPI_MAX, world); + + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + int n = extra_nlen[m]; + double dXInfLocalExtra = 0.0; + for (int i = 0; i < n; i++) + dXInfLocalExtra = MAX (dXInfLocalExtra, fabs (xatom[i])); + double dXInfExtra; + MPI_Allreduce (&dXInfLocalExtra, &dXInfExtra, + 1, MPI_DOUBLE, MPI_MAX, world); + dXInf = MAX (dXInf, dXInfExtra); + } + } + + return( dXInf ); +} + + +/* ---------------------------------------------------------------------- + Private method calc_dot_prod_using_mpi_ +------------------------------------------------------------------------- */ + +double MinHFTN::calc_dot_prod_using_mpi_(const int nIx1, + const int nIx2) const +{ + double dDotLocal = 0.0; + for (int i = 0; i < n3; i++) + dDotLocal += _daAVectors[nIx1][i] * _daAVectors[nIx2][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * i1Atom = _daExtraAtom[nIx1][m]; + double * i2Atom = _daExtraAtom[nIx2][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + dDotLocal += i1Atom[i] * i2Atom[i]; + } + } + + double dDot; + MPI_Allreduce (&dDotLocal, &dDot, 1, MPI_DOUBLE, MPI_SUM, world); + + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) { + double * i1Global = _daExtraGlobal[nIx1]; + double * i2Global = _daExtraGlobal[nIx2]; + dDot += i1Global[i] * i2Global[i]; + } + } + + return( dDot ); +} + +/* ---------------------------------------------------------------------- + Private method calc_grad_dot_v_using_mpi_ +------------------------------------------------------------------------- */ + +double MinHFTN::calc_grad_dot_v_using_mpi_(const int nIx) const +{ + //---- ASSUME THAT FORCES HAVE BEEN EVALUATED AT DESIRED ATOM POSITIONS. + //---- REMEMBER THAT FORCES = -GRADIENT. + + double dGradDotVLocal = 0.0; + for (int i = 0; i < n3; i++) + dGradDotVLocal += - _daAVectors[nIx][i] * atom->f[0][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * iAtom = _daExtraAtom[nIx][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + dGradDotVLocal += - iAtom[i] * fextra_atom[m][i]; + } + } + + double dGradDotV; + MPI_Allreduce (&dGradDotVLocal, &dGradDotV, 1, MPI_DOUBLE, MPI_SUM, world); + + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) { + double * iGlobal = _daExtraGlobal[nIx]; + dGradDotV += - iGlobal[i] * fextra[i]; + } + } + + return( dGradDotV ); +} + +/* ---------------------------------------------------------------------- + Private method calc_dhd_dd_using_mpi_ +------------------------------------------------------------------------- */ + +void MinHFTN::calc_dhd_dd_using_mpi_(double & dDHD, + double & dDD) const +{ + double dDHDLocal = 0.0; + double dDDLocal = 0.0; + for (int i = 0; i < n3; i++) { + dDHDLocal += _daAVectors[VEC_CG_D][i] * _daAVectors[VEC_CG_HD][i]; + dDDLocal += _daAVectors[VEC_CG_D][i] * _daAVectors[VEC_CG_D][i]; + } + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * dAtom = _daExtraAtom[VEC_CG_D][m]; + double * hdAtom = _daExtraAtom[VEC_CG_HD][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) { + dDHDLocal += dAtom[i] * hdAtom[i]; + dDDLocal += dAtom[i] * dAtom[i]; + } + } + } + + double daDotsLocal[2]; + daDotsLocal[0] = dDHDLocal; + daDotsLocal[1] = dDDLocal; + double daDots[2]; + MPI_Allreduce (daDotsLocal, daDots, 2, MPI_DOUBLE, MPI_SUM, world); + + if (nextra_global) { + double * dGlobal = _daExtraGlobal[VEC_CG_D]; + double * hdGlobal = _daExtraGlobal[VEC_CG_HD]; + for (int i = 0; i < nextra_global; i++) { + daDots[0] += dGlobal[i] * hdGlobal[i]; + daDots[1] += dGlobal[i] * dGlobal[i]; + } + } + dDHD = daDots[0]; + dDD = daDots[1]; + + return; +} + +/* ---------------------------------------------------------------------- + Private method calc_ppnew_pdold_using_mpi_ +------------------------------------------------------------------------- */ + +void MinHFTN::calc_ppnew_pdold_using_mpi_(double & dPnewDotPnew, + double & dPoldDotD) const +{ + double dPnewDotPnewLocal = 0.0; + double dPoldDotDLocal = 0.0; + for (int i = 0; i < n3; i++) { + dPnewDotPnewLocal + += _daAVectors[VEC_CG_P][i] * _daAVectors[VEC_CG_P][i]; + dPoldDotDLocal + += _daAVectors[VEC_DIF1][i] * _daAVectors[VEC_CG_D][i]; + } + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * dAtom = _daExtraAtom[VEC_CG_D][m]; + double * pAtom = _daExtraAtom[VEC_CG_P][m]; + double * d1Atom = _daExtraAtom[VEC_DIF1][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) { + dPnewDotPnewLocal += pAtom[i] * pAtom[i]; + dPoldDotDLocal += d1Atom[i] * dAtom[i]; + } + } + } + + double daDotsLocal[2]; + daDotsLocal[0] = dPnewDotPnewLocal; + daDotsLocal[1] = dPoldDotDLocal; + double daDots[2]; + MPI_Allreduce (daDotsLocal, daDots, 2, MPI_DOUBLE, MPI_SUM, world); + + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) { + double * dGlobal = _daExtraGlobal[VEC_CG_D]; + double * pGlobal = _daExtraGlobal[VEC_CG_P]; + double * d1Global = _daExtraGlobal[VEC_DIF1]; + daDots[0] += pGlobal[i] * pGlobal[i]; + daDots[1] += d1Global[i] * dGlobal[i]; + } + } + + dPnewDotPnew = daDots[0]; + dPoldDotD = daDots[1]; + + return; +} + +/* ---------------------------------------------------------------------- + Private method calc_plengths_using_mpi_ +------------------------------------------------------------------------- */ + +void MinHFTN::calc_plengths_using_mpi_(double & dStepLength2, + double & dStepLengthInf) const +{ + double dPPLocal = 0.0; + double dPInfLocal = 0.0; + for (int i = 0; i < n3; i++) { + dPPLocal += _daAVectors[VEC_CG_P][i] * _daAVectors[VEC_CG_P][i]; + dPInfLocal = MAX (dPInfLocal, fabs (_daAVectors[VEC_CG_P][i])); + } + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * pAtom = _daExtraAtom[VEC_CG_P][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) { + dPPLocal += pAtom[i] * pAtom[i]; + dPInfLocal = MAX (dPInfLocal, fabs (pAtom[i])); + } + } + } + + double dPP; + MPI_Allreduce (&dPPLocal, &dPP, 1, MPI_DOUBLE, MPI_SUM, world); + + double dPInf; + MPI_Allreduce (&dPInfLocal, &dPInf, 1, MPI_DOUBLE, MPI_MAX, world); + + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) { + double * pGlobal = _daExtraGlobal[VEC_CG_P]; + dPP += pGlobal[i] * pGlobal[i]; + dPInf = MAX (dPInf, fabs (pGlobal[i])); + } + } + + dStepLength2 = sqrt (dPP); + dStepLengthInf = dPInf; + return; +} + +/* ---------------------------------------------------------------------- + Private method step_exceeds_TR_ +------------------------------------------------------------------------- */ + +bool MinHFTN::step_exceeds_TR_(const double dTrustRadius, + const double dPP, + const double dPD, + const double dDD, + double & dTau) const +{ + double dPnewNorm2; + double dPnewNormInf; + calc_plengths_using_mpi_ (dPnewNorm2, dPnewNormInf); + + if (dPnewNorm2 > dTrustRadius) { + dTau = compute_to_tr_ (dPP, dPD, dDD, dTrustRadius, + false, 0.0, 0.0, 0.0); + return( true ); + } + + //---- STEP LENGTH IS NOT TOO LONG. + dTau = 0.0; + return( false ); +} + +/* ---------------------------------------------------------------------- + Private method step_exceeds_DMAX_ + + Check that atoms do not move too far: + for atom coordinates: limit is dmax + for extra per-atom DOF: limit is extra_max[] + for extra global DOF: limit is set by modify->max_alpha, + which calls fix_box_relax->max_alpha +------------------------------------------------------------------------- */ + +bool MinHFTN::step_exceeds_DMAX_(void) const +{ + double dAlpha = dmax * sqrt (_nNumUnknowns); + + double dPInfLocal = 0.0; + for (int i = 0; i < n3; i++) + dPInfLocal = MAX (dPInfLocal, fabs (_daAVectors[VEC_CG_P][i])); + double dPInf; + MPI_Allreduce (&dPInfLocal, &dPInf, 1, MPI_DOUBLE, MPI_MAX, world); + if (dPInf > dmax) + return( true ); + if (dPInf > MACHINE_EPS) + dAlpha = MIN (dAlpha, dmax / dPInf); + + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * pAtom = _daExtraAtom[VEC_CG_P][m]; + dPInfLocal = 0.0; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + dPInfLocal = MAX (dPInfLocal, fabs (pAtom[i])); + MPI_Allreduce (&dPInfLocal, &dPInf, 1, MPI_DOUBLE, MPI_MAX, world); + if (dPInf > extra_max[m]) + return( true ); + if (dPInf > MACHINE_EPS) + dAlpha = MIN (dAlpha, extra_max[m] / dPInf); + } + } + + if (nextra_global) { + //---- IF THE MAXIMUM DISTANCE THAT THE GLOBAL BOX CONSTRAINT WILL + //---- ALLOW IS SMALLER THAN THE PROPOSED DISTANCE, THEN THE STEP + //---- IS TOO LONG. PROPOSED DISTANCE IS ESTIMATED BY |P|_INF. + double dAlphaExtra = modify->max_alpha (_daExtraGlobal[VEC_CG_P]); + if (dAlphaExtra < dAlpha) + return( true ); + } + + //---- STEP LENGTH IS NOT TOO LONG. + return( false ); +} + +/* ---------------------------------------------------------------------- + Private method adjust_step_to_tau_ + Adjust the step so that VEC_CG_P = VEC_DIF1 + tau * VEC_CG_D. +------------------------------------------------------------------------- */ + +void MinHFTN::adjust_step_to_tau_(const double tau) +{ + if (nextra_global) { + double * pGlobal = _daExtraGlobal[VEC_CG_P]; + double * dGlobal = _daExtraGlobal[VEC_CG_D]; + double * d1Global = _daExtraGlobal[VEC_DIF1]; + for (int i = 0; i < nextra_global; i++) + pGlobal[i] = d1Global[i] + (tau * dGlobal[i]); + } + for (int i = 0; i < n3; i++) { + _daAVectors[VEC_CG_P][i] = _daAVectors[VEC_DIF1][i] + + (tau * _daAVectors[VEC_CG_D][i]); + } + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * pAtom = _daExtraAtom[VEC_CG_P][m]; + double * dAtom = _daExtraAtom[VEC_CG_D][m]; + double * d1Atom = _daExtraAtom[VEC_DIF1][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + pAtom[i] = d1Atom[i] + (tau * dAtom[i]); + } + } + return; +} + +/* ---------------------------------------------------------------------- + Private method compute_to_tr_ + Return the value tau that solves + || p + tau*d ||_2 = dTrustRadius + + If both roots are considered, the TR method chooses the one that minimizes + grad^T (p + tau*d) + 0.5 (p + tau*d)^T H (p + tau*d) + + @param[in] dPP - p^T p + @param[in] dPD - p^T d + @param[in] dDD - d^T d + @param[in] dTrustRadius - distance to match + @param[in] bConsiderBothRoots - evaluate both roots, or return the positive + @param[in] dDHD - d^T H d + @param[in] dPdotHD - p^T H d + @param[in] dGradDotD - grad(x_k)^T d +------------------------------------------------------------------------- */ + +double MinHFTN::compute_to_tr_(const double dPP, + const double dPD, + const double dDD, + const double dTrustRadius, + const bool bConsiderBothRoots, + const double dDHD, + const double dPdotHD, + const double dGradDotD) const +{ + //---- SOLVE A QUADRATIC EQUATION FOR TAU. + //---- THE COEFFICIENTS ARE SUCH THAT THERE ARE ALWAYS TWO REAL ROOTS, + //---- ONE POSITIVE AND ONE NEGATIVE. + + //---- CHECK FOR ERRONEOUS DATA. + if ( (dDD <= 0.0) || (dPP < 0.0) || (dTrustRadius < 0.0) + || (dTrustRadius * dTrustRadius < dPP) ) { + printf ("HFTN internal error - bad data given to compute_to_tr_()\n"); + return( 0.0 ); + } + + double dTRsqrd = dTrustRadius * dTrustRadius; + double dDiscr = (dPD * dPD) - (dDD * (dPP - dTRsqrd)); + dDiscr = MAX (0.0, dDiscr); //-- SHOULD NEVER BE NEGATIVE + dDiscr = sqrt (dDiscr); + + double dRootPos = (-dPD + dDiscr) / dDD; + double dRootNeg = (-dPD - dDiscr) / dDD; + + if (bConsiderBothRoots == false) + return( dRootPos ); + + //---- EVALUATE THE CG OBJECTIVE FUNCTION FOR EACH ROOT. + double dTmpTerm = dGradDotD + dPdotHD; + double dCgRedPos = (dRootPos * dTmpTerm) + (0.5 * dRootPos*dRootPos * dDHD); + double dCgRedNeg = (dRootNeg * dTmpTerm) + (0.5 * dRootNeg*dRootNeg * dDHD); + + if ((-dCgRedPos) > (-dCgRedNeg)) + return( dRootPos ); + else + return( dRootNeg ); +} + +/* ---------------------------------------------------------------------- + Private method evaluate_dir_der_ + Compute the directional derivative using a finite difference approximation. + This is equivalent to the Hessian times direction vector p. + As a side effect, the method computes the energy and forces at x. + + On input these must be defined: + atom->x - positions at x + atom->f - ignored + nIxDir - VEC_ index of the direction p + nIxResult - ignored + On output these are defined: + atom->x - unchanged + atom->f - forces evaluated at x, only if bEvaluateAtX is true + nIxDir - unchanged + nIxResult - directional derivative Hp + During processing these are modified: + VEC_DIF1 + VEC_DIF2 + + @param[in] bUseForwardDiffs - if true use forward difference approximation, + else use central difference + @param[in] nIxDir - VEC_ index of the direction + @param[in] nIxResult - VEC_ index to place the result + (it is acceptable for nIxDir = nIxResult) + @param[in] bEvaluateAtX - if true, then evaluate at x before returning + @param[out] dNewEnergy - energy at x, if bEvaluateAtX is true + @param[out] dNewForce2 - |F|_2 at x, if bEvaluateAtX is true +------------------------------------------------------------------------- */ + +void MinHFTN::evaluate_dir_der_(const bool bUseForwardDiffs, + const int nIxDir, + const int nIxResult, + const bool bEvaluateAtX, + double & dNewEnergy) +{ + //---- COMPUTE THE MAGNITUDE OF THE DIRECTION VECTOR: |p|_2. + double dDirNorm2SqrdLocal = 0.0; + for (int i = 0; i < n3; i++) + dDirNorm2SqrdLocal + += _daAVectors[nIxDir][i] * _daAVectors[nIxDir][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * iAtom = _daExtraAtom[nIxDir][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + dDirNorm2SqrdLocal += iAtom[i] * iAtom[i]; + } + } + double dDirNorm2Sqrd = 0.0; + MPI_Allreduce (&dDirNorm2SqrdLocal, &dDirNorm2Sqrd, + 1, MPI_DOUBLE, MPI_SUM, world); + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) { + double * iGlobal = _daExtraGlobal[nIxDir]; + dDirNorm2Sqrd += iGlobal[i] * iGlobal[i]; + } + } + double dDirNorm2 = sqrt (dDirNorm2Sqrd); + + //---- IF THE STEP IS TOO SMALL, RETURN ZERO FOR THE DERIVATIVE. + if (dDirNorm2 == 0.0) { + for (int i = 0; i < n3; i++) + _daAVectors[nIxResult][i] = 0.0; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * iAtom = _daExtraAtom[nIxDir][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + iAtom[i] = 0; + } + } + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) + _daExtraGlobal[nIxDir][i] = 0.0; + } + + if (bEvaluateAtX) { + dNewEnergy = energy_force (0); + neval++; + } + + return; + } + + //---- FORWARD DIFFERENCES ARE LESS ACCURATE THAN CENTRAL DIFFERENCES, + //---- BUT REQUIRE ONLY 2 ENERGY+FORCE EVALUATIONS VERSUS 3 EVALUATIONS. + //---- STORAGE REQUIREMENTS ARE THE SAME. + + if (bUseForwardDiffs) { + //---- EQUATION IS FROM THE OLD LAMMPS VERSION, SAND98-8201. + double dEps = 2.0 * sqrt (1000.0 * MACHINE_EPS) / dDirNorm2; + + //---- SAVE A COPY OF x. + fix_minimize->store_box(); + for (int i = 0; i < n3; i++) + _daAVectors[VEC_DIF1][i] = atom->x[0][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * d1Atom = _daExtraAtom[VEC_DIF1][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + d1Atom[i] = xatom[i]; + } + } + if (nextra_global) { + modify->min_pushstore(); + modify->min_store(); + } + + //---- EVALUATE FORCES AT x + eps*p. + if (nextra_global) + modify->min_step (dEps, _daExtraGlobal[nIxDir]); + for (int i = 0; i < n3; i++) + atom->x[0][i] += dEps * _daAVectors[nIxDir][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * iAtom = _daExtraAtom[nIxDir][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xatom[i] += dEps * iAtom[i]; + } + } + energy_force (0); + neval++; + + //---- STORE THE FORCE IN DIF2. + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) + _daExtraGlobal[VEC_DIF2][i] = fextra[i]; + } + for (int i = 0; i < n3; i++) + _daAVectors[VEC_DIF2][i] = atom->f[0][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * d2Atom = _daExtraAtom[VEC_DIF2][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + d2Atom[i] = fextra_atom[m][i]; + } + } + + //---- MOVE BACK TO x AND EVALUATE FORCES. + if (nextra_global) { + modify->min_step (0.0, _daExtraGlobal[VEC_DIF1]); + modify->min_popstore(); + } + for (int i = 0; i < n3; i++) + atom->x[0][i] = _daAVectors[VEC_DIF1][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * d1Atom = _daExtraAtom[VEC_DIF1][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xatom[i] += d1Atom[i]; + } + } + dNewEnergy = energy_force (0); + neval++; + + //---- COMPUTE THE DIFFERENCE VECTOR: [grad(x + eps*p) - grad(x)] / eps. + //---- REMEMBER THAT FORCES = -GRADIENT. + for (int i = 0; i < n3; i++) + _daAVectors[nIxResult][i] + = (atom->f[0][i] - _daAVectors[VEC_DIF2][i]) / dEps; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * iAtom = _daExtraAtom[nIxResult][m]; + double * d2Atom = _daExtraAtom[VEC_DIF2][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + iAtom[i] = (fextra_atom[m][i] - d2Atom[i]) / dEps; + } + } + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) + _daExtraGlobal[nIxResult][i] + = (fextra[i] - _daExtraGlobal[VEC_DIF2][i]) / dEps; + } + } + + else { //-- bUseForwardDiffs == false + //---- EQUATION IS FROM THE OLD LAMMPS VERSION, SAND98-8201. + double dEps = pow (3000.0 * MACHINE_EPS, 0.33333333) / dDirNorm2; + + //---- SAVE A COPY OF x. + fix_minimize->store_box(); + for (int i = 0; i < n3; i++) + _daAVectors[VEC_DIF1][i] = atom->x[0][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * d1Atom = _daExtraAtom[VEC_DIF1][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + d1Atom[i] = xatom[i]; + } + } + if (nextra_global) { + modify->min_pushstore(); + modify->min_store(); + } + + //---- EVALUATE FORCES AT x + eps*p. + if (nextra_global) + modify->min_step (dEps, _daExtraGlobal[nIxDir]); + for (int i = 0; i < n3; i++) + atom->x[0][i] += dEps * _daAVectors[nIxDir][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * iAtom = _daExtraAtom[nIxDir][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xatom[i] += dEps * iAtom[i]; + } + } + energy_force (0); + neval++; + + //---- STORE THE FORCE IN DIF2. + if (nextra_global) { + for (int i = 0; i < nextra_global; i++) + _daExtraGlobal[VEC_DIF2][i] = fextra[i]; + } + for (int i = 0; i < n3; i++) + _daAVectors[VEC_DIF2][i] = atom->f[0][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * d2Atom = _daExtraAtom[VEC_DIF2][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + d2Atom[i] = fextra_atom[m][i]; + } + } + + //---- EVALUATE FORCES AT x - eps*p. + if (nextra_global) + modify->min_step (-dEps, _daExtraGlobal[nIxDir]); + for (int i = 0; i < n3; i++) + atom->x[0][i] = _daAVectors[VEC_DIF1][i] + - dEps * _daAVectors[nIxDir][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * iAtom = _daExtraAtom[nIxDir][m]; + double * d1Atom = _daExtraAtom[VEC_DIF1][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xatom[i] = d1Atom[i] - dEps * iAtom[i]; + } + } + energy_force (0); + neval++; + + //---- COMPUTE THE DIFFERENCE VECTOR: + //---- [grad(x + eps*p) - grad(x - eps*p)] / 2*eps. + //---- REMEMBER THAT FORCES = -GRADIENT. + if (nextra_global) { + double * iGlobal = _daExtraGlobal[nIxResult]; + double * d2Global = _daExtraGlobal[VEC_DIF2]; + for (int i = 0; i < nextra_global; i++) + iGlobal[i] = (fextra[i] - d2Global[i]) / (2.0 + dEps); + } + for (int i = 0; i < n3; i++) + _daAVectors[nIxResult][i] + = (atom->f[0][i] - _daAVectors[VEC_DIF2][i]) / (2.0 * dEps); + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * iAtom = _daExtraAtom[nIxResult][m]; + double * d2Atom = _daExtraAtom[VEC_DIF2][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + iAtom[i] = (fextra_atom[m][i] - d2Atom[i]) / (2.0 + dEps); + } + } + + if (bEvaluateAtX) { + //---- EVALUATE FORCES AT x. + if (nextra_global) { + modify->min_step (0.0, _daExtraGlobal[VEC_DIF1]); + modify->min_popstore(); + } + for (int i = 0; i < n3; i++) + atom->x[0][i] = _daAVectors[VEC_DIF1][i]; + if (nextra_atom) { + for (int m = 0; m < nextra_atom; m++) { + double * xatom = xextra_atom[m]; + double * d1Atom = _daExtraAtom[VEC_DIF1][m]; + int n = extra_nlen[m]; + for (int i = 0; i < n; i++) + xatom[i] = d1Atom[i]; + } + } + dNewEnergy = energy_force (0); + neval++; + } + } + + return; +} + +/* ---------------------------------------------------------------------- + Private method open_hftn_print_file_ +------------------------------------------------------------------------- */ + +void MinHFTN::open_hftn_print_file_(void) +{ + int nMyRank; + MPI_Comm_rank (world, &nMyRank); + + char szTmp[50]; + sprintf (szTmp, "progress_MinHFTN_%d.txt", nMyRank); + _fpPrint = fopen (szTmp, "w"); + if (_fpPrint == NULL) { + printf ("*** MinHFTN cannot open file '%s'\n", szTmp); + printf ("*** continuing...\n"); + return; + } + + fprintf (_fpPrint, " Iter Evals Energy |F|_2" + " Step TR used |step|_2 ared pred\n"); + return; +} + +/* ---------------------------------------------------------------------- + Private method hftn_print_line_ + Step types: + 1 - Nw (inner iteration converged like a Newton step) + 2 - TR (inner iteration reached the trust region boundary) + 3 - Neg (inner iteration ended with negative curvature) +------------------------------------------------------------------------- */ + +void MinHFTN::hftn_print_line_(const bool bIsStepAccepted, + const int nIteration, + const int nTotalEvals, + const double dEnergy, + const double dForce2, + const int nStepType, + const double dTrustRadius, + const double dStepLength2, + const double dActualRed, + const double dPredictedRed) const +{ + const char sFormat1[] + = " %4d %5d %14.8f %11.5e\n"; + const char sFormatA[] + = " %4d %5d %14.8f %11.5e %3s %9.3e %8.2e %10.3e %10.3e\n"; + const char sFormatR[] + = "r %4d %5d %14.8f %11.5e %3s %9.3e %8.2e %10.3e %10.3e\n"; + + if (_fpPrint == NULL) + return; + + char sStepType[4]; + if (nStepType == NO_CGSTEP_BECAUSE_F_TOL_SATISFIED) + strcpy (sStepType, " - "); + else if (nStepType == CGSTEP_NEWTON) + strcpy (sStepType, "Nw "); + else if (nStepType == CGSTEP_TO_TR) + strcpy (sStepType, "TR "); + else if (nStepType == CGSTEP_TO_DMAX) + strcpy (sStepType, "dmx"); + else if (nStepType == CGSTEP_NEGATIVE_CURVATURE) + strcpy (sStepType, "Neg"); + else if (nStepType == CGSTEP_MAX_INNER_ITERS) + strcpy (sStepType, "its"); + else + strcpy (sStepType, "???"); + + if (nIteration == -1) { + fprintf (_fpPrint, sFormat1, + 0, nTotalEvals, dEnergy, dForce2); + } + else { + if (bIsStepAccepted) + fprintf (_fpPrint, sFormatA, + nIteration, nTotalEvals, dEnergy, dForce2, + sStepType, dTrustRadius, dStepLength2, + dActualRed, dPredictedRed); + else + fprintf (_fpPrint, sFormatR, + nIteration, nTotalEvals, dEnergy, dForce2, + sStepType, dTrustRadius, dStepLength2, + dActualRed, dPredictedRed); + } + + fflush (_fpPrint); + return; +} + +/* ---------------------------------------------------------------------- + Private method close_hftn_print_file_ +------------------------------------------------------------------------- */ + +void MinHFTN::close_hftn_print_file_(void) +{ + if (_fpPrint != NULL) fclose (_fpPrint); + return; +} diff -Naur lammps-28Sep09/src/min_hftn.h lammps-29Sep09/src/min_hftn.h --- lammps-28Sep09/src/min_hftn.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-29Sep09/src/min_hftn.h 2009-09-29 08:10:33.000000000 -0600 @@ -0,0 +1,119 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef MIN_HFTN_H +#define MIN_HFTN_H + +#include "min.h" + +namespace LAMMPS_NS +{ + //---- THE ALGORITHM NEEDS TO STORE THIS MANY ATOM-BASED VECTORS, + //---- IN ADDITION TO ATOM POSITIONS AND THE FORCE VECTOR. + + static const int NUM_HFTN_ATOM_BASED_VECTORS = 7; + static const int VEC_XK = 0; //-- ATOM POSITIONS AT SUBITER START + static const int VEC_CG_P = 1; //-- STEP p IN CG SUBITER + static const int VEC_CG_D = 2; //-- DIRECTION d IN CG SUBITER + static const int VEC_CG_HD = 3; //-- HESSIAN-VECTOR PRODUCT Hd + static const int VEC_CG_R = 4; //-- RESIDUAL r IN CG SUBITER + static const int VEC_DIF1 = 5; //-- FOR FINITE DIFFERENCING + static const int VEC_DIF2 = 6; //-- FOR FINITE DIFFERENCING + +class MinHFTN : public Min +{ + public: + + MinHFTN (LAMMPS *); + ~MinHFTN (void); + + void init_style(); + void setup_style(); + void reset_vectors(); + int iterate (int); + + private: + + //---- ATOM-BASED STORAGE VECTORS. + double * _daAVectors[NUM_HFTN_ATOM_BASED_VECTORS]; + double ** _daExtraAtom[NUM_HFTN_ATOM_BASED_VECTORS]; + + //---- GLOBAL DOF STORAGE. ELEMENT [0] IS X0 (XK), NOT USED IN THIS ARRAY. + double * _daExtraGlobal[NUM_HFTN_ATOM_BASED_VECTORS]; + + int _nNumUnknowns; + FILE * _fpPrint; + + int execute_hftn_ (const bool bPrintProgress, + const double dInitialEnergy, + const double dInitialForce2, + double & dFinalEnergy, + double & dFinalForce2); + bool compute_inner_cg_step_ (const double dTrustRadius, + const double dForceTol, + const int nMaxEvals, + const bool bHaveEvalAtXin, + const double dEnergyAtXin, + const double dForce2AtXin, + double & dEnergyAtXout, + double & dForce2AtXout, + int & nStepType, + double & dStepLength2, + double & dStepLengthInf); + double calc_xinf_using_mpi_ (void) const; + double calc_dot_prod_using_mpi_ (const int nIx1, + const int nIx2) const; + double calc_grad_dot_v_using_mpi_ (const int nIx) const; + void calc_dhd_dd_using_mpi_ (double & dDHD, + double & dDD) const; + void calc_ppnew_pdold_using_mpi_ (double & dPnewDotPnew, + double & dPoldDotD) const; + void calc_plengths_using_mpi_ (double & dStepLength2, + double & dStepLengthInf) const; + bool step_exceeds_TR_ (const double dTrustRadius, + const double dPP, + const double dPD, + const double dDD, + double & dTau) const; + bool step_exceeds_DMAX_ (void) const; + void adjust_step_to_tau_ (const double tau); + double compute_to_tr_ (const double dPP, + const double dPD, + const double dDD, + const double dTrustRadius, + const bool bConsiderBothRoots, + const double dDHD, + const double dPdotHD, + const double dGradDotD) const; + void evaluate_dir_der_ (const bool bUseForwardDiffs, + const int nIxDir, + const int nIxResult, + const bool bEvaluateAtX, + double & dNewEnergy); + void open_hftn_print_file_ (void); + void hftn_print_line_ (const bool bIsStepAccepted, + const int nIteration, + const int nTotalEvals, + const double dEnergy, + const double dForce2, + const int nStepType, + const double dTrustRadius, + const double dStepLength2, + const double dActualRed, + const double dPredictedRed) const; + void close_hftn_print_file_ (void); +}; + +} + +#endif diff -Naur lammps-28Sep09/src/modify.cpp lammps-29Sep09/src/modify.cpp --- lammps-28Sep09/src/modify.cpp 2009-08-28 14:28:23.000000000 -0600 +++ lammps-29Sep09/src/modify.cpp 2009-09-29 08:10:33.000000000 -0600 @@ -441,6 +441,28 @@ } /* ---------------------------------------------------------------------- + mange state of extra dof on a stack, only for relevant fixes +------------------------------------------------------------------------- */ + +void Modify::min_clearstore() +{ + for (int i = 0; i < n_min_energy; i++) + fix[list_min_energy[i]]->min_clearstore(); +} + +void Modify::min_pushstore() +{ + for (int i = 0; i < n_min_energy; i++) + fix[list_min_energy[i]]->min_pushstore(); +} + +void Modify::min_popstore() +{ + for (int i = 0; i < n_min_energy; i++) + fix[list_min_energy[i]]->min_popstore(); +} + +/* ---------------------------------------------------------------------- displace extra dof along vector hextra, only for relevant fixes ------------------------------------------------------------------------- */ diff -Naur lammps-28Sep09/src/modify.h lammps-29Sep09/src/modify.h --- lammps-28Sep09/src/modify.h 2009-08-28 14:28:23.000000000 -0600 +++ lammps-29Sep09/src/modify.h 2009-09-29 08:10:33.000000000 -0600 @@ -65,6 +65,9 @@ double min_energy(double *); void min_store(); void min_step(double, double *); + void min_clearstore(); + void min_pushstore(); + void min_popstore(); double max_alpha(double *); int min_dof(); diff -Naur lammps-28Sep09/src/style.h lammps-29Sep09/src/style.h --- lammps-28Sep09/src/style.h 2009-09-25 09:15:47.000000000 -0600 +++ lammps-29Sep09/src/style.h 2009-09-29 08:10:33.000000000 -0600 @@ -297,11 +297,13 @@ #ifdef MinimizeInclude #include "min_cg.h" +#include "min_hftn.h" #include "min_sd.h" #endif #ifdef MinimizeClass MinimizeStyle(cg,MinCG) +MinimizeStyle(hftn,MinHFTN) MinimizeStyle(sd,MinSD) # endif