ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13713 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2015-07-24 17:57:41 +00:00
parent 052f644d13
commit c6c44b5fe0
24 changed files with 2221 additions and 73 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
doc/Manual.html
View File

@ -199,13 +199,13 @@ it gives quick access to documentation for all LAMMPS commands.
<BR> <BR>
6.21 <A HREF = "Section_howto.html#howto_21">Calculating viscosity</A> 6.21 <A HREF = "Section_howto.html#howto_21">Calculating viscosity</A>
<BR> <BR>
6.22 <A HREF = "howto_22">Calculating a diffusion coefficient</A> 6.22 <A HREF = "Section_howto.html#howto_22">Calculating a diffusion coefficient</A>
<BR> <BR>
6.23 <A HREF = "howto_23">Using chunks to calculate system properties</A> 6.23 <A HREF = "Section_howto.html#howto_23">Using chunks to calculate system properties</A>
<BR> <BR>
6.24 <A HREF = "howto_24">Setting parameters for pppm/disp</A> 6.24 <A HREF = "Section_howto.html#howto_24">Setting parameters for pppm/disp</A>
<BR> <BR>
6.25 <A HREF = "howto_25">Adiabatic core/shell model</A> 6.25 <A HREF = "Section_howto.html#howto_25">Adiabatic core/shell model</A>
<BR></UL> <BR></UL>
<LI><A HREF = "Section_example.html">Example problems</A> <LI><A HREF = "Section_example.html">Example problems</A>
@ -412,6 +412,16 @@ it gives quick access to documentation for all LAMMPS commands.

5
doc/Manual.txt
View File

@ -228,6 +228,11 @@ it gives quick access to documentation for all LAMMPS commands.
:link(howto_19,Section_howto.html#howto_19) :link(howto_19,Section_howto.html#howto_19)
:link(howto_20,Section_howto.html#howto_20) :link(howto_20,Section_howto.html#howto_20)
:link(howto_21,Section_howto.html#howto_21) :link(howto_21,Section_howto.html#howto_21)
:link(howto_22,Section_howto.html#howto_22)
:link(howto_23,Section_howto.html#howto_23)
:link(howto_24,Section_howto.html#howto_24)
:link(howto_25,Section_howto.html#howto_25)
:link(howto_26,Section_howto.html#howto_26)
:link(mod_1,Section_modify.html#mod_1) :link(mod_1,Section_modify.html#mod_1)
:link(mod_2,Section_modify.html#mod_2) :link(mod_2,Section_modify.html#mod_2)

69
doc/Section_commands.html
View File

@ -433,12 +433,13 @@ g = GPU, i = USER-INTEL, k = KOKKOS, o = USER-OMP, t = OPT.
package</A>. package</A>.
</P> </P>
<DIV ALIGN=center><TABLE BORDER=1 > <DIV ALIGN=center><TABLE BORDER=1 >
<TR ALIGN="center"><TD ><A HREF = "fix_adapt_fep.html">adapt/fep</A></TD><TD ><A HREF = "fix_addtorque.html">addtorque</A></TD><TD ><A HREF = "fix_atc.html">atc</A></TD><TD ><A HREF = "fix_ave_spatial_sphere.html">ave/spatial/sphere</A></TD><TD ><A HREF = "fix_colvars.html">colvars</A></TD><TD ><A HREF = "fix_gle.html">gle</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "fix_adapt_fep.html">adapt/fep</A></TD><TD ><A HREF = "fix_addtorque.html">addtorque</A></TD><TD ><A HREF = "fix_atc.html">atc</A></TD><TD ><A HREF = "fix_ave_spatial_sphere.html">ave/spatial/sphere</A></TD><TD ><A HREF = "fix_drude.html">drude</A></TD><TD ><A HREF = "fix_drude_transform.html">drude/transform/direct</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_imd.html">imd</A></TD><TD ><A HREF = "fix_ipi.html">ipi</A></TD><TD ><A HREF = "fix_langevin_eff.html">langevin/eff</A></TD><TD ><A HREF = "fix_lb_fluid.html">lb/fluid</A></TD><TD ><A HREF = "fix_lb_momentum.html">lb/momentum</A></TD><TD ><A HREF = "fix_lb_pc.html">lb/pc</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "fix_drude_transform.html">drude/transform/reverse</A></TD><TD ><A HREF = "fix_colvars.html">colvars</A></TD><TD ><A HREF = "fix_gle.html">gle</A></TD><TD ><A HREF = "fix_imd.html">imd</A></TD><TD ><A HREF = "fix_ipi.html">ipi</A></TD><TD ><A HREF = "fix_langevin_drude.html">langevin/drude</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_lb_rigid_pc_sphere.html">lb/rigid/pc/sphere</A></TD><TD ><A HREF = "fix_lb_viscous.html">lb/viscous</A></TD><TD ><A HREF = "fix_meso.html">meso</A></TD><TD ><A HREF = "fix_meso_stationary.html">meso/stationary</A></TD><TD ><A HREF = "fix_nh_eff.html">nph/eff</A></TD><TD ><A HREF = "fix_nh_eff.html">npt/eff</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "fix_langevin_eff.html">langevin/eff</A></TD><TD ><A HREF = "fix_lb_fluid.html">lb/fluid</A></TD><TD ><A HREF = "fix_lb_momentum.html">lb/momentum</A></TD><TD ><A HREF = "fix_lb_pc.html">lb/pc</A></TD><TD ><A HREF = "fix_lb_rigid_pc_sphere.html">lb/rigid/pc/sphere</A></TD><TD ><A HREF = "fix_lb_viscous.html">lb/viscous</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_nve_eff.html">nve/eff</A></TD><TD ><A HREF = "fix_nh_eff.html">nvt/eff</A></TD><TD ><A HREF = "fix_nvt_sllod_eff.html">nvt/sllod/eff</A></TD><TD ><A HREF = "fix_phonon.html">phonon</A></TD><TD ><A HREF = "fix_pimd.html">pimd</A></TD><TD ><A HREF = "fix_qbmsst.html">qbmsst</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "fix_meso.html">meso</A></TD><TD ><A HREF = "fix_meso_stationary.html">meso/stationary</A></TD><TD ><A HREF = "fix_nh_eff.html">nph/eff</A></TD><TD ><A HREF = "fix_nh_eff.html">npt/eff</A></TD><TD ><A HREF = "fix_nve_eff.html">nve/eff</A></TD><TD ><A HREF = "fix_nh_eff.html">nvt/eff</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_qeq_reax.html">qeq/reax</A></TD><TD ><A HREF = "fix_qmmm.html">qmmm</A></TD><TD ><A HREF = "fix_qtb.html">qtb</A></TD><TD ><A HREF = "fix_reax_bonds.html">reax/c/bonds</A></TD><TD ><A HREF = "fix_reaxc_species.html">reax/c/species</A></TD><TD ><A HREF = "fix_saed_vtk.html">saed/vtk</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "fix_nvt_sllod_eff.html">nvt/sllod/eff</A></TD><TD ><A HREF = "fix_phonon.html">phonon</A></TD><TD ><A HREF = "fix_pimd.html">pimd</A></TD><TD ><A HREF = "fix_qbmsst.html">qbmsst</A></TD><TD ><A HREF = "fix_qeq_reax.html">qeq/reax</A></TD><TD ><A HREF = "fix_qmmm.html">qmmm</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_smd.html">smd</A></TD><TD ><A HREF = "fix_temp_rescale_eff.html">temp/rescale/eff</A></TD><TD ><A HREF = "fix_ti_rs.html">ti/rs</A></TD><TD ><A HREF = "fix_ti_spring.html">ti/spring</A></TD><TD ><A HREF = "fix_ttm.html">ttm/mod</A> <TR ALIGN="center"><TD ><A HREF = "fix_qtb.html">qtb</A></TD><TD ><A HREF = "fix_reax_bonds.html">reax/c/bonds</A></TD><TD ><A HREF = "fix_reaxc_species.html">reax/c/species</A></TD><TD ><A HREF = "fix_saed_vtk.html">saed/vtk</A></TD><TD ><A HREF = "fix_smd.html">smd</A></TD><TD ><A HREF = "fix_temp_rescale_eff.html">temp/rescale/eff</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_ti_rs.html">ti/rs</A></TD><TD ><A HREF = "fix_ti_spring.html">ti/spring</A></TD><TD ><A HREF = "fix_ttm.html">ttm/mod</A>
</TD></TR></TABLE></DIV> </TD></TR></TABLE></DIV>
<HR> <HR>
@ -473,8 +474,8 @@ package</A>.
</P> </P>
<DIV ALIGN=center><TABLE BORDER=1 > <DIV ALIGN=center><TABLE BORDER=1 >
<TR ALIGN="center"><TD ><A HREF = "compute_ackland_atom.html">ackland/atom</A></TD><TD ><A HREF = "compute_basal_atom.html">basal/atom</A></TD><TD ><A HREF = "compute_fep.html">fep</A></TD><TD ><A HREF = "compute_ke_eff.html">ke/eff</A></TD><TD ><A HREF = "compute_ke_atom_eff.html">ke/atom/eff</A></TD><TD ><A HREF = "compute_meso_e_atom.html">meso_e/atom</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "compute_ackland_atom.html">ackland/atom</A></TD><TD ><A HREF = "compute_basal_atom.html">basal/atom</A></TD><TD ><A HREF = "compute_fep.html">fep</A></TD><TD ><A HREF = "compute_ke_eff.html">ke/eff</A></TD><TD ><A HREF = "compute_ke_atom_eff.html">ke/atom/eff</A></TD><TD ><A HREF = "compute_meso_e_atom.html">meso_e/atom</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "compute_meso_rho_atom.html">meso_rho/atom</A></TD><TD ><A HREF = "compute_meso_t_atom.html">meso_t/atom</A></TD><TD ><A HREF = "compute_saed.html">saed</A></TD><TD ><A HREF = "compute_temp_eff.html">temp/eff</A></TD><TD ><A HREF = "compute_temp_deform_eff.html">temp/deform/eff</A></TD><TD ><A HREF = "compute_temp_region_eff.html">temp/region/eff</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "compute_meso_rho_atom.html">meso_rho/atom</A></TD><TD ><A HREF = "compute_meso_t_atom.html">meso_t/atom</A></TD><TD ><A HREF = "compute_saed.html">saed</A></TD><TD ><A HREF = "compute_temp_drude.html">temp/drude</A></TD><TD ><A HREF = "compute_temp_eff.html">temp/eff</A></TD><TD ><A HREF = "compute_temp_deform_eff.html">temp/deform/eff</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "compute_temp_rotate.html">temp/rotate</A></TD><TD ><A HREF = "compute_xrd.html">xrd</A> <TR ALIGN="center"><TD ><A HREF = "compute_temp_region_eff.html">temp/region/eff</A></TD><TD ><A HREF = "compute_temp_rotate.html">temp/rotate</A></TD><TD ><A HREF = "compute_xrd.html">xrd</A>
</TD></TR></TABLE></DIV> </TD></TR></TABLE></DIV>
<HR> <HR>
@ -492,30 +493,31 @@ KOKKOS, o = USER-OMP, t = OPT.
<DIV ALIGN=center><TABLE BORDER=1 > <DIV ALIGN=center><TABLE BORDER=1 >
<TR ALIGN="center"><TD ><A HREF = "pair_none.html">none</A></TD><TD ><A HREF = "pair_hybrid.html">hybrid</A></TD><TD ><A HREF = "pair_hybrid.html">hybrid/overlay</A></TD><TD ><A HREF = "pair_adp.html">adp (o)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_none.html">none</A></TD><TD ><A HREF = "pair_hybrid.html">hybrid</A></TD><TD ><A HREF = "pair_hybrid.html">hybrid/overlay</A></TD><TD ><A HREF = "pair_adp.html">adp (o)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_airebo.html">airebo (o)</A></TD><TD ><A HREF = "pair_beck.html">beck (go)</A></TD><TD ><A HREF = "pair_body.html">body</A></TD><TD ><A HREF = "pair_bop.html">bop</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_airebo.html">airebo (o)</A></TD><TD ><A HREF = "pair_beck.html">beck (go)</A></TD><TD ><A HREF = "pair_body.html">body</A></TD><TD ><A HREF = "pair_bop.html">bop</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_born.html">born (go)</A></TD><TD ><A HREF = "pair_born.html">born/coul/long (cgo)</A></TD><TD ><A HREF = "pair_born.html">born/coul/msm (o)</A></TD><TD ><A HREF = "pair_born.html">born/coul/wolf (go)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_born.html">born (go)</A></TD><TD ><A HREF = "pair_born.html">born/coul/long (cgo)</A></TD><TD ><A HREF = "pair_born.html">born/coul/long/cs</A></TD><TD ><A HREF = "pair_born.html">born/coul/msm (o)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_brownian.html">brownian (o)</A></TD><TD ><A HREF = "pair_brownian.html">brownian/poly (o)</A></TD><TD ><A HREF = "pair_buck.html">buck (cgko)</A></TD><TD ><A HREF = "pair_buck.html">buck/coul/cut (cgko)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_born.html">born/coul/wolf (go)</A></TD><TD ><A HREF = "pair_brownian.html">brownian (o)</A></TD><TD ><A HREF = "pair_brownian.html">brownian/poly (o)</A></TD><TD ><A HREF = "pair_buck.html">buck (cgko)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_buck.html">buck/coul/long (cgko)</A></TD><TD ><A HREF = "pair_buck.html">buck/coul/msm (o)</A></TD><TD ><A HREF = "pair_buck_long.html">buck/long/coul/long (o)</A></TD><TD ><A HREF = "pair_colloid.html">colloid (go)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_buck.html">buck/coul/cut (cgko)</A></TD><TD ><A HREF = "pair_buck.html">buck/coul/long (cgko)</A></TD><TD ><A HREF = "pair_buck.html">buck/coul/long/cs</A></TD><TD ><A HREF = "pair_buck.html">buck/coul/msm (o)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_comb.html">comb (o)</A></TD><TD ><A HREF = "pair_comb.html">comb3</A></TD><TD ><A HREF = "pair_coul.html">coul/cut (gko)</A></TD><TD ><A HREF = "pair_coul.html">coul/debye (gko)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_buck_long.html">buck/long/coul/long (o)</A></TD><TD ><A HREF = "pair_colloid.html">colloid (go)</A></TD><TD ><A HREF = "pair_comb.html">comb (o)</A></TD><TD ><A HREF = "pair_comb.html">comb3</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_coul.html">coul/dsf (gko)</A></TD><TD ><A HREF = "pair_coul.html">coul/long (gko)</A></TD><TD ><A HREF = "pair_coul.html">coul/msm</A></TD><TD ><A HREF = "pair_coul.html">coul/streitz</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_coul.html">coul/cut (gko)</A></TD><TD ><A HREF = "pair_coul.html">coul/debye (gko)</A></TD><TD ><A HREF = "pair_coul.html">coul/dsf (gko)</A></TD><TD ><A HREF = "pair_coul.html">coul/long (gko)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_coul.html">coul/wolf (ko)</A></TD><TD ><A HREF = "pair_dpd.html">dpd (o)</A></TD><TD ><A HREF = "pair_dpd.html">dpd/tstat (o)</A></TD><TD ><A HREF = "pair_dsmc.html">dsmc</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_coul.html">coul/long/cs</A></TD><TD ><A HREF = "pair_coul.html">coul/msm</A></TD><TD ><A HREF = "pair_coul.html">coul/streitz</A></TD><TD ><A HREF = "pair_coul.html">coul/wolf (ko)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_eam.html">eam (cgkot)</A></TD><TD ><A HREF = "pair_eam.html">eam/alloy (cgkot)</A></TD><TD ><A HREF = "pair_eam.html">eam/fs (cgkot)</A></TD><TD ><A HREF = "pair_eim.html">eim (o)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_dpd.html">dpd (o)</A></TD><TD ><A HREF = "pair_dpd.html">dpd/tstat (o)</A></TD><TD ><A HREF = "pair_dsmc.html">dsmc</A></TD><TD ><A HREF = "pair_eam.html">eam (cgkot)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_gauss.html">gauss (go)</A></TD><TD ><A HREF = "pair_gayberne.html">gayberne (gio)</A></TD><TD ><A HREF = "pair_gran.html">gran/hertz/history (o)</A></TD><TD ><A HREF = "pair_gran.html">gran/hooke (co)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_eam.html">eam/alloy (cgkot)</A></TD><TD ><A HREF = "pair_eam.html">eam/fs (cgkot)</A></TD><TD ><A HREF = "pair_eim.html">eim (o)</A></TD><TD ><A HREF = "pair_gauss.html">gauss (go)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_gran.html">gran/hooke/history (o)</A></TD><TD ><A HREF = "pair_hbond_dreiding.html">hbond/dreiding/lj (o)</A></TD><TD ><A HREF = "pair_hbond_dreiding.html">hbond/dreiding/morse (o)</A></TD><TD ><A HREF = "pair_kim.html">kim</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_gayberne.html">gayberne (gio)</A></TD><TD ><A HREF = "pair_gran.html">gran/hertz/history (o)</A></TD><TD ><A HREF = "pair_gran.html">gran/hooke (co)</A></TD><TD ><A HREF = "pair_gran.html">gran/hooke/history (o)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_lcbop.html">lcbop</A></TD><TD ><A HREF = "pair_line_lj.html">line/lj (o)</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm (cko)</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm/implicit (cko)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_hbond_dreiding.html">hbond/dreiding/lj (o)</A></TD><TD ><A HREF = "pair_hbond_dreiding.html">hbond/dreiding/morse (o)</A></TD><TD ><A HREF = "pair_kim.html">kim</A></TD><TD ><A HREF = "pair_lcbop.html">lcbop</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/long (cgiko)</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/msm</A></TD><TD ><A HREF = "pair_class2.html">lj/class2 (cgko)</A></TD><TD ><A HREF = "pair_class2.html">lj/class2/coul/cut (cko)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_line_lj.html">line/lj (o)</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm (cko)</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/charmm/implicit (cko)</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/long (cgiko)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_class2.html">lj/class2/coul/long (cgko)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut (cgikot)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/cut (cgko)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/debye (cgko)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/msm</A></TD><TD ><A HREF = "pair_class2.html">lj/class2 (cgko)</A></TD><TD ><A HREF = "pair_class2.html">lj/class2/coul/cut (cko)</A></TD><TD ><A HREF = "pair_class2.html">lj/class2/coul/long (cgko)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut/coul/dsf (gko)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/long (cgikot)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/msm (go)</A></TD><TD ><A HREF = "pair_dipole.html">lj/cut/dipole/cut (go)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut (cgikot)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/cut (cgko)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/debye (cgko)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/dsf (gko)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_dipole.html">lj/cut/dipole/long</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/tip4p/cut (o)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/tip4p/long (ot)</A></TD><TD ><A HREF = "pair_lj_expand.html">lj/expand (cgko)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut/coul/long (cgikot)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/msm (go)</A></TD><TD ><A HREF = "pair_dipole.html">lj/cut/dipole/cut (go)</A></TD><TD ><A HREF = "pair_dipole.html">lj/cut/dipole/long</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_gromacs.html">lj/gromacs (cgko)</A></TD><TD ><A HREF = "pair_gromacs.html">lj/gromacs/coul/gromacs (cko)</A></TD><TD ><A HREF = "pair_lj_long.html">lj/long/coul/long (o)</A></TD><TD ><A HREF = "pair_dipole.html">lj/long/dipole/long</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut/tip4p/cut (o)</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/tip4p/long (ot)</A></TD><TD ><A HREF = "pair_lj_expand.html">lj/expand (cgko)</A></TD><TD ><A HREF = "pair_gromacs.html">lj/gromacs (cgko)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_lj_long.html">lj/long/tip4p/long</A></TD><TD ><A HREF = "pair_lj_smooth.html">lj/smooth (co)</A></TD><TD ><A HREF = "pair_lj_smooth_linear.html">lj/smooth/linear (o)</A></TD><TD ><A HREF = "pair_lj96.html">lj96/cut (cgo)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_gromacs.html">lj/gromacs/coul/gromacs (cko)</A></TD><TD ><A HREF = "pair_lj_long.html">lj/long/coul/long (o)</A></TD><TD ><A HREF = "pair_dipole.html">lj/long/dipole/long</A></TD><TD ><A HREF = "pair_lj_long.html">lj/long/tip4p/long</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_lubricate.html">lubricate (o)</A></TD><TD ><A HREF = "pair_lubricate.html">lubricate/poly (o)</A></TD><TD ><A HREF = "pair_lubricateU.html">lubricateU</A></TD><TD ><A HREF = "pair_lubricateU.html">lubricateU/poly</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_lj_smooth.html">lj/smooth (co)</A></TD><TD ><A HREF = "pair_lj_smooth_linear.html">lj/smooth/linear (o)</A></TD><TD ><A HREF = "pair_lj96.html">lj96/cut (cgo)</A></TD><TD ><A HREF = "pair_lubricate.html">lubricate (o)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_meam.html">meam (o)</A></TD><TD ><A HREF = "pair_mie.html">mie/cut (o)</A></TD><TD ><A HREF = "pair_morse.html">morse (cgot)</A></TD><TD ><A HREF = "pair_nb3b_harmonic.html">nb3b/harmonic (o)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_lubricate.html">lubricate/poly (o)</A></TD><TD ><A HREF = "pair_lubricateU.html">lubricateU</A></TD><TD ><A HREF = "pair_lubricateU.html">lubricateU/poly</A></TD><TD ><A HREF = "pair_meam.html">meam (o)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_nm.html">nm/cut (o)</A></TD><TD ><A HREF = "pair_nm.html">nm/cut/coul/cut (o)</A></TD><TD ><A HREF = "pair_nm.html">nm/cut/coul/long (o)</A></TD><TD ><A HREF = "pair_peri.html">peri/eps</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_mie.html">mie/cut (o)</A></TD><TD ><A HREF = "pair_morse.html">morse (cgot)</A></TD><TD ><A HREF = "pair_nb3b_harmonic.html">nb3b/harmonic (o)</A></TD><TD ><A HREF = "pair_nm.html">nm/cut (o)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_peri.html">peri/lps (o)</A></TD><TD ><A HREF = "pair_peri.html">peri/pmb (o)</A></TD><TD ><A HREF = "pair_peri.html">peri/ves</A></TD><TD ><A HREF = "pair_reax.html">reax</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_nm.html">nm/cut/coul/cut (o)</A></TD><TD ><A HREF = "pair_nm.html">nm/cut/coul/long (o)</A></TD><TD ><A HREF = "pair_peri.html">peri/eps</A></TD><TD ><A HREF = "pair_peri.html">peri/lps (o)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_airebo.html">rebo (o)</A></TD><TD ><A HREF = "pair_resquared.html">resquared (go)</A></TD><TD ><A HREF = "pair_snap.html">snap</A></TD><TD ><A HREF = "pair_soft.html">soft (go)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_peri.html">peri/pmb (o)</A></TD><TD ><A HREF = "pair_peri.html">peri/ves</A></TD><TD ><A HREF = "pair_reax.html">reax</A></TD><TD ><A HREF = "pair_airebo.html">rebo (o)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_sw.html">sw (cgkio)</A></TD><TD ><A HREF = "pair_table.html">table (gko)</A></TD><TD ><A HREF = "pair_tersoff.html">tersoff (cko)</A></TD><TD ><A HREF = "pair_tersoff_mod.html">tersoff/mod (ko)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_resquared.html">resquared (go)</A></TD><TD ><A HREF = "pair_snap.html">snap</A></TD><TD ><A HREF = "pair_soft.html">soft (go)</A></TD><TD ><A HREF = "pair_sw.html">sw (cgkio)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_tersoff_zbl.html">tersoff/zbl (ko)</A></TD><TD ><A HREF = "pair_coul.html">tip4p/cut (o)</A></TD><TD ><A HREF = "pair_coul.html">tip4p/long (o)</A></TD><TD ><A HREF = "pair_tri_lj.html">tri/lj (o)</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_table.html">table (gko)</A></TD><TD ><A HREF = "pair_tersoff.html">tersoff (cko)</A></TD><TD ><A HREF = "pair_tersoff_mod.html">tersoff/mod (ko)</A></TD><TD ><A HREF = "pair_tersoff_zbl.html">tersoff/zbl (ko)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_yukawa.html">yukawa (go)</A></TD><TD ><A HREF = "pair_yukawa_colloid.html">yukawa/colloid (go)</A></TD><TD ><A HREF = "pair_zbl.html">zbl (o)</A> <TR ALIGN="center"><TD ><A HREF = "pair_coul.html">tip4p/cut (o)</A></TD><TD ><A HREF = "pair_coul.html">tip4p/long (o)</A></TD><TD ><A HREF = "pair_tri_lj.html">tri/lj (o)</A></TD><TD ><A HREF = "pair_yukawa.html">yukawa (go)</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_yukawa_colloid.html">yukawa/colloid (go)</A></TD><TD ><A HREF = "pair_zbl.html">zbl (o)</A>
</TD></TR></TABLE></DIV> </TD></TR></TABLE></DIV>
<P>These are additional pair styles in USER packages, which can be used <P>These are additional pair styles in USER packages, which can be used
@ -530,7 +532,8 @@ package</A>.
<TR ALIGN="center"><TD ><A HREF = "pair_sdk.html">lj/sdk/coul/long (go)</A></TD><TD ><A HREF = "pair_sdk.html">lj/sdk/coul/msm (o)</A></TD><TD ><A HREF = "pair_lj_sf.html">lj/sf (o)</A></TD><TD ><A HREF = "pair_meam_spline.html">meam/spline</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_sdk.html">lj/sdk/coul/long (go)</A></TD><TD ><A HREF = "pair_sdk.html">lj/sdk/coul/msm (o)</A></TD><TD ><A HREF = "pair_lj_sf.html">lj/sf (o)</A></TD><TD ><A HREF = "pair_meam_spline.html">meam/spline</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_meam_sw_spline.html">meam/sw/spline</A></TD><TD ><A HREF = "pair_quip.html">quip</A></TD><TD ><A HREF = "pair_reax_c.html">reax/c</A></TD><TD ><A HREF = "pair_sph_heatconduction.html">sph/heatconduction</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_meam_sw_spline.html">meam/sw/spline</A></TD><TD ><A HREF = "pair_quip.html">quip</A></TD><TD ><A HREF = "pair_reax_c.html">reax/c</A></TD><TD ><A HREF = "pair_sph_heatconduction.html">sph/heatconduction</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_sph_idealgas.html">sph/idealgas</A></TD><TD ><A HREF = "pair_sph_lj.html">sph/lj</A></TD><TD ><A HREF = "pair_sph_rhosum.html">sph/rhosum</A></TD><TD ><A HREF = "pair_sph_taitwater.html">sph/taitwater</A></TD></TR> <TR ALIGN="center"><TD ><A HREF = "pair_sph_idealgas.html">sph/idealgas</A></TD><TD ><A HREF = "pair_sph_lj.html">sph/lj</A></TD><TD ><A HREF = "pair_sph_rhosum.html">sph/rhosum</A></TD><TD ><A HREF = "pair_sph_taitwater.html">sph/taitwater</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_sph_taitwater_morris.html">sph/taitwater/morris</A></TD><TD ><A HREF = "pair_srp.html">srp</A></TD><TD ><A HREF = "pair_tersoff.html">tersoff/table (o)</A></TD><TD ><A HREF = "pair_lj_soft.html">tip4p/long/soft (o)</A> <TR ALIGN="center"><TD ><A HREF = "pair_sph_taitwater_morris.html">sph/taitwater/morris</A></TD><TD ><A HREF = "pair_srp.html">srp</A></TD><TD ><A HREF = "pair_tersoff.html">tersoff/table (o)</A></TD><TD ><A HREF = "pair_thole.html">thole</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "pair_lj_soft.html">tip4p/long/soft (o)</A>
</TD></TR></TABLE></DIV> </TD></TR></TABLE></DIV>
<HR> <HR>

9
doc/Section_commands.txt
View File

@ -602,10 +602,14 @@ package"_Section_start.html#start_3.
"addtorque"_fix_addtorque.html, "addtorque"_fix_addtorque.html,
"atc"_fix_atc.html, "atc"_fix_atc.html,
"ave/spatial/sphere"_fix_ave_spatial_sphere.html, "ave/spatial/sphere"_fix_ave_spatial_sphere.html,
"drude"_fix_drude.html,
"drude/transform/direct"_fix_drude_transform.html,
"drude/transform/reverse"_fix_drude_transform.html,
"colvars"_fix_colvars.html, "colvars"_fix_colvars.html,
"gle"_fix_gle.html, "gle"_fix_gle.html,
"imd"_fix_imd.html, "imd"_fix_imd.html,
"ipi"_fix_ipi.html, "ipi"_fix_ipi.html,
"langevin/drude"_fix_langevin_drude.html,
"langevin/eff"_fix_langevin_eff.html, "langevin/eff"_fix_langevin_eff.html,
"lb/fluid"_fix_lb_fluid.html, "lb/fluid"_fix_lb_fluid.html,
"lb/momentum"_fix_lb_momentum.html, "lb/momentum"_fix_lb_momentum.html,
@ -726,6 +730,7 @@ package"_Section_start.html#start_3.
"meso_rho/atom"_compute_meso_rho_atom.html, "meso_rho/atom"_compute_meso_rho_atom.html,
"meso_t/atom"_compute_meso_t_atom.html, "meso_t/atom"_compute_meso_t_atom.html,
"saed"_compute_saed.html, "saed"_compute_saed.html,
"temp/drude"_compute_temp_drude.html,
"temp/eff"_compute_temp_eff.html, "temp/eff"_compute_temp_eff.html,
"temp/deform/eff"_compute_temp_deform_eff.html, "temp/deform/eff"_compute_temp_deform_eff.html,
"temp/region/eff"_compute_temp_region_eff.html, "temp/region/eff"_compute_temp_region_eff.html,
@ -754,6 +759,7 @@ KOKKOS, o = USER-OMP, t = OPT.
"bop"_pair_bop.html, "bop"_pair_bop.html,
"born (go)"_pair_born.html, "born (go)"_pair_born.html,
"born/coul/long (cgo)"_pair_born.html, "born/coul/long (cgo)"_pair_born.html,
"born/coul/long/cs"_pair_born.html,
"born/coul/msm (o)"_pair_born.html, "born/coul/msm (o)"_pair_born.html,
"born/coul/wolf (go)"_pair_born.html, "born/coul/wolf (go)"_pair_born.html,
"brownian (o)"_pair_brownian.html, "brownian (o)"_pair_brownian.html,
@ -761,6 +767,7 @@ KOKKOS, o = USER-OMP, t = OPT.
"buck (cgko)"_pair_buck.html, "buck (cgko)"_pair_buck.html,
"buck/coul/cut (cgko)"_pair_buck.html, "buck/coul/cut (cgko)"_pair_buck.html,
"buck/coul/long (cgko)"_pair_buck.html, "buck/coul/long (cgko)"_pair_buck.html,
"buck/coul/long/cs"_pair_buck.html,
"buck/coul/msm (o)"_pair_buck.html, "buck/coul/msm (o)"_pair_buck.html,
"buck/long/coul/long (o)"_pair_buck_long.html, "buck/long/coul/long (o)"_pair_buck_long.html,
"colloid (go)"_pair_colloid.html, "colloid (go)"_pair_colloid.html,
@ -770,6 +777,7 @@ KOKKOS, o = USER-OMP, t = OPT.
"coul/debye (gko)"_pair_coul.html, "coul/debye (gko)"_pair_coul.html,
"coul/dsf (gko)"_pair_coul.html, "coul/dsf (gko)"_pair_coul.html,
"coul/long (gko)"_pair_coul.html, "coul/long (gko)"_pair_coul.html,
"coul/long/cs"_pair_coul.html,
"coul/msm"_pair_coul.html, "coul/msm"_pair_coul.html,
"coul/streitz"_pair_coul.html, "coul/streitz"_pair_coul.html,
"coul/wolf (ko)"_pair_coul.html, "coul/wolf (ko)"_pair_coul.html,
@ -883,6 +891,7 @@ package"_Section_start.html#start_3.
"sph/taitwater/morris"_pair_sph_taitwater_morris.html, "sph/taitwater/morris"_pair_sph_taitwater_morris.html,
"srp"_pair_srp.html, "srp"_pair_srp.html,
"tersoff/table (o)"_pair_tersoff.html, "tersoff/table (o)"_pair_tersoff.html,
"thole"_pair_thole.html,
"tip4p/long/soft (o)"_pair_lj_soft.html :tb(c=4,ea=c) "tip4p/long/soft (o)"_pair_lj_soft.html :tb(c=4,ea=c)
:line :line

9
doc/compute_reduce.html
View File

@ -25,7 +25,7 @@
<I>reduce/region</I> arg = region-ID <I>reduce/region</I> arg = region-ID
region-ID = ID of region to use for choosing atoms region-ID = ID of region to use for choosing atoms
</PRE> </PRE>
<LI>mode = <I>sum</I> or <I>min</I> or <I>max</I> or <I>ave</I> <LI>mode = <I>sum</I> or <I>min</I> or <I>max</I> or <I>ave</I> or <I>sumsq</I> or <I>avesq</I>
<LI>one or more inputs can be listed <LI>one or more inputs can be listed
@ -71,7 +71,12 @@ vectors.
option adds the values in the vector into a global total. The <I>min</I> option adds the values in the vector into a global total. The <I>min</I>
or <I>max</I> options find the minimum or maximum value across all vector or <I>max</I> options find the minimum or maximum value across all vector
values. The <I>ave</I> setting adds the vector values into a global total, values. The <I>ave</I> setting adds the vector values into a global total,
then divides by the number of values in the vector. then divides by the number of values in the vector. The <I>sumsq</I>
option sums the square of the values in the vector into a global
total. The <I>avesq</I> setting does the same as <I>sumsq</I>, then divdes the
sum of squares by the number of values. The last two options can be
useful for calculating the variance of some quantity, e.g. variance =
sumsq - ave^2.
</P> </P>
<P>Each listed input is operated on independently. For per-atom inputs, <P>Each listed input is operated on independently. For per-atom inputs,
the group specified with this command means only atoms within the the group specified with this command means only atoms within the

9
doc/compute_reduce.txt
View File

@ -18,7 +18,7 @@ style = {reduce} or {reduce/region} :l
{reduce} arg = none {reduce} arg = none
{reduce/region} arg = region-ID {reduce/region} arg = region-ID
region-ID = ID of region to use for choosing atoms :pre region-ID = ID of region to use for choosing atoms :pre
mode = {sum} or {min} or {max} or {ave} :l mode = {sum} or {min} or {max} or {ave} or {sumsq} or {avesq} :l
one or more inputs can be listed :l one or more inputs can be listed :l
input = x, y, z, vx, vy, vz, fx, fy, fz, c_ID, c_ID\[N\], f_ID, f_ID\[N\], v_name :l input = x, y, z, vx, vy, vz, fx, fy, fz, c_ID, c_ID\[N\], f_ID, f_ID\[N\], v_name :l
x,y,z,vx,vy,vz,fx,fy,fz = atom attribute (position, velocity, force component) x,y,z,vx,vy,vz,fx,fy,fz = atom attribute (position, velocity, force component)
@ -58,7 +58,12 @@ The reduction operation is specified by the {mode} setting. The {sum}
option adds the values in the vector into a global total. The {min} option adds the values in the vector into a global total. The {min}
or {max} options find the minimum or maximum value across all vector or {max} options find the minimum or maximum value across all vector
values. The {ave} setting adds the vector values into a global total, values. The {ave} setting adds the vector values into a global total,
then divides by the number of values in the vector. then divides by the number of values in the vector. The {sumsq}
option sums the square of the values in the vector into a global
total. The {avesq} setting does the same as {sumsq}, then divdes the
sum of squares by the number of values. The last two options can be
useful for calculating the variance of some quantity, e.g. variance =
sumsq - ave^2.
Each listed input is operated on independently. For per-atom inputs, Each listed input is operated on independently. For per-atom inputs,
the group specified with this command means only atoms within the the group specified with this command means only atoms within the

62
doc/compute_temp_drude.html Normal file
View File

@ -0,0 +1,62 @@
<HTML>
<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A>
</CENTER>
<HR>
<H3>compute temp/drude command
</H3>
<P><B>Syntax:</B>
</P>
<PRE>compute ID group-ID temp/drude
</PRE>
<UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command
<LI>temp/drude = style name of this compute command
</UL>
<P><B>Examples:</B>
</P>
<PRE>compute TDRUDE all temp/drude
</PRE>
<P><B>Description:</B>
</P>
<P>Define a computation that calculates the temperature based on the
center-of-mass velocities of pairs of Drude cores and Drude particles,
bonded by springs. This compute is designed to be used with the
<A HREF = "tutorial_drude.html">thermalized Drude oscillator model</A>. Polarizable
models in LAMMPS are described in <A HREF = "Section_howto.html#howto_25">this
Section</A>.
</P>
<P>Specifically, this compute enables calculation of the temperature of
the Drude particles in relative coordinates with respect to their
cores.
</P>
<P><B>Output info:</B>
</P>
<P>This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1-6.
These values can be used by any command that uses global scalar or
vector values from a compute as input. See <A HREF = "Section_howto.html#howto_15">this
section</A> for an overview of LAMMPS output
options.
</P>
<P>The scalar value calculated by this compute is "intensive". The
vector are "extensive".
</P>
<P>The scalar value will be in temperature <A HREF = "units.html">units</A>. The
vector values will be in energy <A HREF = "units.html">units</A>.
</P>
<P><B>Restrictions:</B>
</P>
<P>The number of core-Drude pairs contributing to the temperature is
assumed to be constant for the duration of the run.
</P>
<P><B>Related commands:</B> none
</P>
<P><B>Default:</B> none
</P>
</HTML>

57
doc/compute_temp_drude.txt Normal file
View File

@ -0,0 +1,57 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
compute temp/drude command :h3
[Syntax:]
compute ID group-ID temp/drude :pre
ID, group-ID are documented in "compute"_compute.html command
temp/drude = style name of this compute command :ul
[Examples:]
compute TDRUDE all temp/drude :pre
[Description:]
Define a computation that calculates the temperature based on the
center-of-mass velocities of pairs of Drude cores and Drude particles,
bonded by springs. This compute is designed to be used with the
"thermalized Drude oscillator model"_tutorial_drude.html. Polarizable
models in LAMMPS are described in "this
Section"_Section_howto.html#howto_25.
Specifically, this compute enables calculation of the temperature of
the Drude particles in relative coordinates with respect to their
cores.
[Output info:]
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1-6.
These values can be used by any command that uses global scalar or
vector values from a compute as input. See "this
section"_Section_howto.html#howto_15 for an overview of LAMMPS output
options.
The scalar value calculated by this compute is "intensive". The
vector are "extensive".
The scalar value will be in temperature "units"_units.html. The
vector values will be in energy "units"_units.html.
[Restrictions:]
The number of core-Drude pairs contributing to the temperature is
assumed to be constant for the duration of the run.
[Related commands:] none
[Default:] none

60
doc/fix_drude.html Normal file
View File

@ -0,0 +1,60 @@
<HTML>
<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A>
</CENTER>
<HR>
<H3>fix drude command
</H3>
<P><B>Syntax:</B>
</P>
<PRE>fix ID group-ID drude flag1 flag2 ... flagN
</PRE>
<UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command
<LI>drude = style name of this fix command
<LI>tag = Drude flag for each atom type (1 to N) in the system
</UL>
<P><B>Examples:</B>
</P>
<PRE>fix 1 all drude 1 1 0 1 0 2 2 2
fix 1 all drude C C N C N D D D
</PRE>
<P><B>Description:</B>
</P>
<P>Assign each atom type in the system to be one of 3 kinds of atoms
within the Drude polarization model. This compute is designed to be
used with the <A HREF = "tutorial_drude.html">thermalized Drude oscillator
model</A>. Polarizable models in LAMMPS
are described in <A HREF = "Section_howto.html#howto_25">this Section</A>.
</P>
<P>The three possible types can be designated with an integer (0,1,2)
or capital letter (N,C,D):
</P>
<UL><LI>0 or N = non-polarizable atom (not part of Drude model)
<LI>1 or C = Drude core
<LI>2 or D = Drude electron
</UL>
<P><B>Restrictions:</B>
</P>
<P>This fix should be invoked before any other commands that implement
the Drude oscillator model, such as <A HREF = "fix_langevin_drude.html">fix
langevin_drude</A>, <A HREF = "fix_drude_transform.html">fix
drude/transform</A>, <A HREF = "compute_temp_drude.html">compute
temp/drude</A>, <A HREF = "pair_thole.html">pair_style
thole</A>.
</P>
<P><B>Related commands:</B>
</P>
<P><A HREF = "fix_langevin_drude.html">fix langevin_drude</A>, <A HREF = "fix_drude_transform.html">fix
drude/transform</A>, <A HREF = "compute_temp_drude.html">compute
temp/drude</A>, <A HREF = "pair_thole.html">pair_style
thole</A>
</P>
<P><B>Default:</B> None
</P>
</HTML>

56
doc/fix_drude.txt Normal file
View File

@ -0,0 +1,56 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix drude command :h3
[Syntax:]
fix ID group-ID drude flag1 flag2 ... flagN :pre
ID, group-ID are documented in "fix"_fix.html command
drude = style name of this fix command
tag = Drude flag for each atom type (1 to N) in the system :ul
[Examples:]
fix 1 all drude 1 1 0 1 0 2 2 2
fix 1 all drude C C N C N D D D :pre
[Description:]
Assign each atom type in the system to be one of 3 kinds of atoms
within the Drude polarization model. This compute is designed to be
used with the "thermalized Drude oscillator
model"_tutorial_drude.html. Polarizable models in LAMMPS
are described in "this Section"_Section_howto.html#howto_25.
The three possible types can be designated with an integer (0,1,2)
or capital letter (N,C,D):
0 or N = non-polarizable atom (not part of Drude model)
1 or C = Drude core
2 or D = Drude electron :ul
[Restrictions:]
This fix should be invoked before any other commands that implement
the Drude oscillator model, such as "fix
langevin_drude"_fix_langevin_drude.html, "fix
drude/transform"_fix_drude_transform.html, "compute
temp/drude"_compute_temp_drude.html, "pair_style
thole"_pair_thole.html.
[Related commands:]
"fix langevin_drude"_fix_langevin_drude.html, "fix
drude/transform"_fix_drude_transform.html, "compute
temp/drude"_compute_temp_drude.html, "pair_style
thole"_pair_thole.html
[Default:] None

195
doc/fix_drude_transform.html Normal file
View File

@ -0,0 +1,195 @@
<HTML>
<script type="text/javascript"
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
</script>
<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A>
</CENTER>
<HR>
<H3>fix drude/transform/direct command
</H3>
<H3>fix drude/transform/inverse command
</H3>
<P><B>Syntax:</B>
</P>
<PRE>fix ID group-ID style keyword value ...
</PRE>
<UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command
<LI>style = <I>drude/transform/direct</I> or <I>drude/transform/inverse</I>
<LI>zero or more keywords may be appended to drude/transform/inverse
<LI>keyword = <I>temp</I>
<PRE> <I>temp</I> value = <I>yes</I> or <I>no</I> = do or do not calculate reduced temperatures of core and Drude particles
</PRE>
</UL>
<P><B>Examples:</B>
</P>
<PRE>fix 3 all drude/transform/direct
fix 1 all drude/transform/inverse temp yes
</PRE>
<P><B>Description:</B>
</P>
<P>Transform the coordinates of Drude oscillators from real to reduced
and back for thermalizing the Drude oscillators as described in
<A HREF = "#Lamoureux">(Lamoureux)</A> using a Nose-Hoover thermostat. This fix is
designed to be used with the <A HREF = "tutorial_drude.html">thermalized Drude oscillator
model</A>. Polarizable models in LAMMPS are
described in <A HREF = "Section_howto.html#howto_25">this Section</A>.
</P>
<P>Drude oscillators are a pair of atoms representing a single
polarizable atom. Ideally, the mass of Drude particles would vanish
and their positions would be determined self-consistently by iterative
minimization of the energy, the cores' positions being fixed. It is
however more efficient and it yields comparable results, if the Drude
oscillators (the motion of the Drude particle relative to the core)
are thermalized at a low temperature. In that case, the Drude
particles need a small mass.
</P>
<P>The thermostats act on the reduced degrees of freedom, which are
defined by the following equations. Note that in these equations
upper case denotes atomic or center of mass values and lower case
denotes Drude particle or dipole values. Primes denote the transformed
(reduced) values, while bare letters denote the original values.
</P>
<P>Masses: \begin{equation} M' = M + m \end{equation}
\begin{equation} m' = \frac {M\, m } {M'} \end{equation}
Positions: \begin{equation} X' = \frac {M\, X + m\, x} {M'}
\end{equation} \begin{equation} x' = x - X \end{equation}
Velocities: \begin{equation} V' = \frac {M\, V + m\, v} {M'}
\end{equation} \begin{equation} v' = v - V \end{equation}
Forces: \begin{equation} F' = F + f \end{equation}
\begin{equation} f' = \frac { M\, f - m\, F} {M'}
\end{equation}
</P>
<P>This transform conserves the total kinetic energy
\begin{equation} \frac 1 2 \, (M\, V^2\ + m\, v^2)
= \frac 1 2 \, (M'\, V'^2\ + m'\, v'^2) \end{equation}
and the virial defined with absolute positions
\begin{equation} X\, F + x\, f = X'\, F' + x'\, f' \end{equation}
</P>
<P>The <I>temp</I> keyword specifies wheterh temperatures in reduced units for
cores and Drude particles are calculated. If the <I>temp</I> option is set
to <I>yes</I> the reduced temperatures, degrees of freedom, and kinetic
energies are calculated and can be accessed as explained below;
otherwise the reduced values are not calculated.
</P>
<HR>
<P>This fix requires each atom know whether it is a Drude particle or
not. You must therefore use the <A HREF = "fix_drude.html">fix drude</A> command to
specify the Drude status of each atom type.
</P>
<P>IMPORTANT NOTE: only the Drude core atoms need to be in the group
specified for this fix. A Drude electron will be transformed together
with its cores even if it is not itself in the group. It is safe to
include Drude electrons or non-polarizable atoms in the group. The
non-polarizable atoms will simply not be transformed.
</P>
<HR>
<P>This fix does NOT perform time integration. It only transform masses,
coordinates, velocities and forces. Thus you must use separate time
integration fixes, like <A HREF = "fix_nve.html">fix nve</A> or <A HREF = "fix_nh.html">fix
npt</A> to actually update the velocities and positions of
atoms. In order to thermalize the reduced degrees of freedom at
different temperatures, two Nose-Hoover thermostats must be defined,
acting on two distinct groups.
</P>
<P>IMPORTANT NOTE: The <I>fix drude/transform/direct</I> command must appear
before any Nose-Hoover thermostating fixes. The <I>fix
drude/transform/inverse</I> command must appear after any Nose-Hoover
thermostating fixes.
</P>
<P>Example:
</P>
<PRE>fix fDIRECT all drude/transform/direct
fix fNVT gCORES nvt temp 300.0 300.0 100.0
fix fNVT gDRUDES nvt temp 1.0 1.0 100.0
fix fINVERSE all drude/transform/inverse
</PRE>
<P>In this example, <I>gCORES</I> is the group of the atom cores and <I>gDRUDES</I>
is the group of the Drude particles (electrons). The centers of mass
of the Drude oscillators will be thermostated at 300.0 and the
internal degrees of freedom will be thermostated at 1.0.
</P>
<P>In addition, if you want to use a barostat to simulate a system at
constant pressure, only one of the Nose-Hoover fixes must be <I>npt</I>,
the other one should be <I>nvt</I>. You must add a <I>compute temp/com</I> and a
<I>fix_modify</I> command so that the temperature of the <I>npt</I> fix be just
that of its group but the pressure be the overall pressure
<I>thermo_press</I>.
</P>
Example:
<BR>
<PRE>compute cTEMP_CORE gCORES temp/com
fix fDIRECT all drude/transform/direct
fix fNPT gCORES npt temp 298.0 298.0 100.0 iso 1.0 1.0 500.0
fix_modify fNPT temp cTEMP_CORE press thermo_press
fix fNVT gDRUDES nvt temp 5.0 5.0 100.0
fix fINVERSE all drude/transform/inverse
</PRE>
<P>In this example, <I>gCORES</I> is the group of the atom cores and <I>gDRUDES</I>
is the group of the Drude particles. The centers of mass of the Drude
oscillators will be thermostated at 298.0 and the internal degrees of
freedom will be thermostated at 5.0. The whole system will be
barostated at 1.0.
</P>
<P>In order to avoid the flying ice cube problem (irreversible transfer
of linear momentum to the center of mass of the system), you may need
to add a <I>fix momentum</I> command like such as
</P>
<PRE>fix fMOMENTUM all momentum 100 linear 1 1 1
</PRE>
<HR>
<P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
</P>
<P>If the <I>temp yes</I> keyword is used in <I>fix drude/transform/inverse</I>
this fix computes a global vector with 6 components which can be
accessed by various <A HREF = "Section_howto.html#howto_15">output commands</A>. The
meaning of the components are
</P>
<OL><LI>temperature of the centers of mass (temperature units)
<LI>temperature of the dipoles (temperature units)
<LI>number of degrees of freedom of the centers of mass
<LI>number of degrees of freedom of the dipoles
<LI>kinetic energy of the centers of mass (energy units)
<LI>kinetic energy of the dipoles (energy units)
</OL>
<P>No information about this fix is written to <A HREF = "restart.html">binary restart
files</A>.
</P>
<P><B>Restrictions:</B> none
</P>
<P><B>Related commands:</B>
</P>
<P><A HREF = "fix_drude.html">fix drude</A>,
<A HREF = "fix_langevin_drude.html">fix langevin/drude</A>,
<A HREF = "compute_temp_drude.html">compute temp/drude</A>,
<A HREF = "pair_thole.html">pair_style thole</A>
</P>
<P><B>Default:</B>
</P>
<P>The option defaults are temp = no.
</P>
<HR>
<A NAME = "Lamoureux"></A>
<P><B>(Lamoureux)</B> Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003).
</P>
</HTML>

183
doc/fix_drude_transform.txt Normal file
View File

@ -0,0 +1,183 @@
<script type="text/javascript"
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
</script>
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix drude/transform/direct command :h3
fix drude/transform/inverse command :h3
[Syntax:]
fix ID group-ID style keyword value ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
style = {drude/transform/direct} or {drude/transform/inverse} :l
zero or more keywords may be appended to drude/transform/inverse :l
keyword = {temp} :l
{temp} value = {yes} or {no} = do or do not calculate reduced temperatures of core and Drude particles :pre
:ule
[Examples:]
fix 3 all drude/transform/direct
fix 1 all drude/transform/inverse temp yes :pre
[Description:]
Transform the coordinates of Drude oscillators from real to reduced
and back for thermalizing the Drude oscillators as described in
"(Lamoureux)"_#Lamoureux using a Nose-Hoover thermostat. This fix is
designed to be used with the "thermalized Drude oscillator
model"_tutorial_drude.html. Polarizable models in LAMMPS are
described in "this Section"_Section_howto.html#howto_25.
Drude oscillators are a pair of atoms representing a single
polarizable atom. Ideally, the mass of Drude particles would vanish
and their positions would be determined self-consistently by iterative
minimization of the energy, the cores' positions being fixed. It is
however more efficient and it yields comparable results, if the Drude
oscillators (the motion of the Drude particle relative to the core)
are thermalized at a low temperature. In that case, the Drude
particles need a small mass.
The thermostats act on the reduced degrees of freedom, which are
defined by the following equations. Note that in these equations
upper case denotes atomic or center of mass values and lower case
denotes Drude particle or dipole values. Primes denote the transformed
(reduced) values, while bare letters denote the original values.
Masses: \begin\{equation\} M' = M + m \end\{equation\}
\begin\{equation\} m' = \frac \{M\, m \} \{M'\} \end\{equation\}
Positions: \begin\{equation\} X' = \frac \{M\, X + m\, x\} \{M'\}
\end\{equation\} \begin\{equation\} x' = x - X \end\{equation\}
Velocities: \begin\{equation\} V' = \frac \{M\, V + m\, v\} \{M'\}
\end\{equation\} \begin\{equation\} v' = v - V \end\{equation\}
Forces: \begin\{equation\} F' = F + f \end\{equation\}
\begin\{equation\} f' = \frac \{ M\, f - m\, F\} \{M'\}
\end\{equation\}
This transform conserves the total kinetic energy
\begin\{equation\} \frac 1 2 \, (M\, V^2\ + m\, v^2)
= \frac 1 2 \, (M'\, V'^2\ + m'\, v'^2) \end\{equation\}
and the virial defined with absolute positions
\begin\{equation\} X\, F + x\, f = X'\, F' + x'\, f' \end\{equation\}
The {temp} keyword specifies wheterh temperatures in reduced units for
cores and Drude particles are calculated. If the {temp} option is set
to {yes} the reduced temperatures, degrees of freedom, and kinetic
energies are calculated and can be accessed as explained below;
otherwise the reduced values are not calculated.
:line
This fix requires each atom know whether it is a Drude particle or
not. You must therefore use the "fix drude"_fix_drude.html command to
specify the Drude status of each atom type.
IMPORTANT NOTE: only the Drude core atoms need to be in the group
specified for this fix. A Drude electron will be transformed together
with its cores even if it is not itself in the group. It is safe to
include Drude electrons or non-polarizable atoms in the group. The
non-polarizable atoms will simply not be transformed.
:line
This fix does NOT perform time integration. It only transform masses,
coordinates, velocities and forces. Thus you must use separate time
integration fixes, like "fix nve"_fix_nve.html or "fix
npt"_fix_nh.html to actually update the velocities and positions of
atoms. In order to thermalize the reduced degrees of freedom at
different temperatures, two Nose-Hoover thermostats must be defined,
acting on two distinct groups.
IMPORTANT NOTE: The {fix drude/transform/direct} command must appear
before any Nose-Hoover thermostating fixes. The {fix
drude/transform/inverse} command must appear after any Nose-Hoover
thermostating fixes.
Example:
fix fDIRECT all drude/transform/direct
fix fNVT gCORES nvt temp 300.0 300.0 100.0
fix fNVT gDRUDES nvt temp 1.0 1.0 100.0
fix fINVERSE all drude/transform/inverse :pre
In this example, {gCORES} is the group of the atom cores and {gDRUDES}
is the group of the Drude particles (electrons). The centers of mass
of the Drude oscillators will be thermostated at 300.0 and the
internal degrees of freedom will be thermostated at 1.0.
In addition, if you want to use a barostat to simulate a system at
constant pressure, only one of the Nose-Hoover fixes must be {npt},
the other one should be {nvt}. You must add a {compute temp/com} and a
{fix_modify} command so that the temperature of the {npt} fix be just
that of its group but the pressure be the overall pressure
{thermo_press}.
Example: :b
compute cTEMP_CORE gCORES temp/com
fix fDIRECT all drude/transform/direct
fix fNPT gCORES npt temp 298.0 298.0 100.0 iso 1.0 1.0 500.0
fix_modify fNPT temp cTEMP_CORE press thermo_press
fix fNVT gDRUDES nvt temp 5.0 5.0 100.0
fix fINVERSE all drude/transform/inverse :pre
In this example, {gCORES} is the group of the atom cores and {gDRUDES}
is the group of the Drude particles. The centers of mass of the Drude
oscillators will be thermostated at 298.0 and the internal degrees of
freedom will be thermostated at 5.0. The whole system will be
barostated at 1.0.
In order to avoid the flying ice cube problem (irreversible transfer
of linear momentum to the center of mass of the system), you may need
to add a {fix momentum} command like such as
fix fMOMENTUM all momentum 100 linear 1 1 1 :pre
:line
[Restart, fix_modify, output, run start/stop, minimize info:]
If the {temp yes} keyword is used in {fix drude/transform/inverse}
this fix computes a global vector with 6 components which can be
accessed by various "output commands"_Section_howto.html#howto_15. The
meaning of the components are
temperature of the centers of mass (temperature units)
temperature of the dipoles (temperature units)
number of degrees of freedom of the centers of mass
number of degrees of freedom of the dipoles
kinetic energy of the centers of mass (energy units)
kinetic energy of the dipoles (energy units) :ol
No information about this fix is written to "binary restart
files"_restart.html.
[Restrictions:] none
[Related commands:]
"fix drude"_fix_drude.html,
"fix langevin/drude"_fix_langevin_drude.html,
"compute temp/drude"_compute_temp_drude.html,
"pair_style thole"_pair_thole.html
[Default:]
The option defaults are temp = no.
:line
:link(Lamoureux)
[(Lamoureux)] Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003).

272
doc/fix_langevin_drude.html Normal file
View File

@ -0,0 +1,272 @@
<HTML>
<script type="text/javascript"
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
</script>
<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A>
</CENTER>
<HR>
<H3>fix langevin/drude command
</H3>
<P><B>Syntax:</B>
</P>
<PRE>fix ID group-ID langevin/drude Tcom damp_com seed_com Tdrude damp_drude seed_drude keyword values ...
</PRE>
<UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command
<LI>langevin/drude = style name of this fix command
<LI>Tcom = desired temperature of the centers of mass (temperature units)
<LI>damp_com = damping parameter for the thermostat on centers of mass (time units)
<LI>seed_com = random number seed to use for white noise of the thermostat on centers of mass (positive integer)
<LI>Tdrude = desired temperature of the Drude oscillators (temperature units)
<LI>damp_drude = damping parameter for the thermostat on Drude oscillators (time units)
<LI>seed_drude = random number seed to use for white noise of the thermostat on Drude oscillators (positive integer)
<LI>zero or more keyword/value pairs may be appended
<LI>keyword = <I>zero</I>
<PRE> <I>zero</I> value = <I>no</I> or <I>yes</I>
<I>no</I> = do not set total random force on centers of mass to zero
<I>yes</I> = set total random force on centers of mass to zero
</PRE>
</UL>
<P><B>Examples:</B>
</P>
<PRE>fix 3 all langevin/drude 300.0 100.0 19377 1.0 20.0 83451
fix 1 all langevin/drude 298.15 100.0 19377 5.0 10.0 83451 zero yes
</PRE>
<P><B>Description:</B>
</P>
<P>Apply two Langevin thermostats as described in <A HREF = "#Jiang">(Jiang)</A> for
thermalizing the reduced degrees of freedom of Drude oscillators.
This link describes how to use the <A HREF = "tutorial_drude.html">thermalized Drude oscillator
model</A> in LAMMPS and polarizable models in LAMMPS
are discussed in <A HREF = "Section_howto.html#howto_25">this Section</A>.
</P>
<P>Drude oscillators are a way to simulate polarizables atoms, by
splitting them into a core and a Drude particle bound by a harmonic
bond. The thermalization works by transforming the particles degrees
of freedom by these equations. In these equations upper case denotes
atomic or center of mass values and lower case denotes Drude particle
or dipole values. Primes denote the transformed (reduced) values,
while bare letters denote the original values.
</P>
Velocities:
\begin{equation} V' = \frac {M\, V + m\, v} {M'} \end{equation}
\begin{equation} v' = v - V \end{equation}
Masses:
\begin{equation} M' = M + m \end{equation}
\begin{equation} m' = \frac {M\, m } {M'} \end{equation}
The Langevin forces are computed as
\begin{equation} F' = - \frac {M'} {\mathtt{damp\_com}}\, V' + F_r' \end{equation}
\begin{equation} f' = - \frac {m'} {\mathtt{damp\_drude}}\, v' + f_r' \end{equation}
\(F_r'\) is a random force proportional to
\(\sqrt { \frac {2\, k_B \mathtt{Tcom}\, m'}
{\mathrm dt\, \mathtt{damp\_com} }
} \).
<BR>
\(f_r'\) is a random force proportional to
\(\sqrt { \frac {2\, k_B \mathtt{Tdrude}\, m'}
{\mathrm dt\, \mathtt{damp\_drude} }
} \).
<BR>
<P>Then the real forces acting on the particles are computed from the inverse
transform:
\begin{equation} F = \frac M {M'}\, F' - f' \end{equation}
\begin{equation} f = \frac m {M'}\, F' + f' \end{equation}
</P>
<P>For Drude pairs (core + electron), the center of mass and the dipole
are thermostated if (and only if) the core atom is in the specified
group.
</P>
<P>Note that the thermostat effect of this fix is applied to only the
translational degrees of freedom of the particles, which is an
important consideration if finite-size particles, which have
rotational degrees of freedom, are being thermostated. The
translational degrees of freedom can also have a bias velocity removed
from them before thermostating takes place; see the description below.
</P>
<P>IMPORTANT NOTE: Like the <A HREF = "fix_langevin.html">fix langevin</A> command,
this fix does NOT perform time integration. It only modifies forces to
effect thermostating. Thus you must use a separate time integration
fix, like <A HREF = "fix_nve.html">fix nve</A> or <A HREF = "fix_nh.html">fix nph</A> to actually
update the velocities and positions of atoms using the modified
forces. Likewise, this fix should not normally be used on atoms that
also have their temperature controlled by another fix - e.g. by <A HREF = "fix_nh.html">fix
nvt</A> or <A HREF = "fix_temp_rescale.html">fix temp/rescale</A> commands.
</P>
<P>See <A HREF = "Section_howto.html#howto_16">this howto section</A> of the manual for
a discussion of different ways to compute temperature and perform
thermostating.
</P>
<HR>
<P>This fix requires each atom know whether it is a Drude particle or
not. You must therefore use the <A HREF = "fix_drude.html">fix drude</A> command to
specify the Drude status of each atom type.
</P>
<P>IMPORTANT NOTE: only the Drude core atoms need to be in the group
specified for this fix. A Drude electron will be transformed together
with its cores even if it is not itself in the group. It is safe to
include Drude electrons or non-polarizable atoms in the group. The
non-polarizable atoms will simply not be thermostatted as if they had
a massless Drude partner (electron).
</P>
<P>IMPORTANT NOTE: Ghost atoms need to know their velocity for this fix
to act correctly. You must use the <A HREF = "comm_modify.html">comm_modify</A>
command to enable this, e.g.
</P>
<PRE>comm_modify vel yes
</PRE>
<HR>
<P><I>Tcom</I> is the target temperature of the centers of mass, which would
be used to thermostate the non-polarizable atoms. <I>Tdrude</I> is the
(normally low) target temperature of the core-Drude particle pairs
(dipoles). <I>Tcom</I> and <I>Tdrude</I> can be specified as an equal-style
<A HREF = "variable.html">variable</A>. If the value is a variable, it should be
specified as v_name, where name is the variable name. In this case,
the variable will be evaluated each timestep, and its value used to
determine the target temperature.
</P>
<P>Equal-style variables can specify formulas with various mathematical
functions, and include <A HREF = "thermo_style.html">thermo_style</A> command
keywords for the simulation box parameters and timestep and elapsed
time. Thus it is easy to specify a time-dependent temperature.
</P>
<P>Like other fixes that perform thermostating, this fix can be used with
<A HREF = "compute.html">compute commands</A> that remove a "bias" from the atom
velocities. E.g. removing the center-of-mass velocity from a group of
atoms. This is not done by default, but only if the
<A HREF = "fix_modify.html">fix_modify</A> command is used to assign a temperature
compute to this fix that includes such a bias term. See the doc pages
for individual <A HREF = "compute.html">compute commands</A> to determine which ones
include a bias. In this case, the thermostat works in the following
manner: bias is removed from each atom, thermostating is performed on
the remaining thermal degrees of freedom, and the bias is added back
in. NOTE: this feature has not been tested.
</P>
<P>Note: The temperature thermostating the core-Drude particle pairs
should be chosen low enough, so as to mimic as closely as possible the
self-consistent minimization. It must however be high enough, so that
the dipoles can follow the local electric field exerted by the
neighbouring atoms. The optimal value probably depends on the
temperature of the centers of mass and on the mass of the Drude
particles.
</P>
<P><I>damp_com</I> is the characteristic time for reaching thermal equilibrium
of the centers of mass. For example, a value of 100.0 means to relax
the temperature of the centers of mass in a timespan of (roughly) 100
time units (tau or fmsec or psec - see the <A HREF = "units.html">units</A>
command). <I>damp_drude</I> is the characteristic time for reaching
thermal equilibrium of the dipoles. It is typically a few timesteps.
</P>
<P>The number <I>seed_com</I> and <I>seed_drude</I> are positive integers. They set
the seeds of the Marsaglia random number generators used for
generating the random forces on centers of mass and on the
dipoles. Each processor uses the input seed to generate its own unique
seed and its own stream of random numbers. Thus the dynamics of the
system will not be identical on two runs on different numbers of
processors.
</P>
<P>The keyword <I>zero</I> can be used to eliminate drift due to the
thermostat on centers of mass. Because the random forces on different
centers of mass are independent, they do not sum exactly to zero. As
a result, this fix applies a small random force to the entire system,
and the momentum of the total center of mass of the system undergoes a
slow random walk. If the keyword <I>zero</I> is set to <I>yes</I>, the total
random force on the centers of mass is set exactly to zero by
subtracting off an equal part of it from each center of mass in the
group. As a result, the total center of mass of a system with zero
initial momentum will not drift over time.
</P>
<HR>
<P>Example for rigid bodies in the NPT ensemble:
</P>
<PRE>comm_modify vel yes
fix TEMP all langevin/drude 300. 100. 1256 1. 20. 13977 zero yes
fix NPH ATOMS rigid/nph/small molecule iso 1. 1. 500.
fix NVE DRUDES nve
thermo_style custom step cpu etotal ke pe ebond ecoul elong press vol temp c_TATOM c_TDRUDE f_TEMP[1] f_TEMP[2]
</PRE>
<P>Comments:
</P>
<UL><LI>Drude particles should not be in the rigid group, otherwise the Drude oscillators will be frozen and the system will lose its polarizability.
<LI><I>zero yes</I> avoids a drift of the center of mass of the system, but is a bit slower.
<LI>use two different random seeds to avoid unphysical correlations.
<LI>temperature is controlled by the fix <I>langevin/drude</I>, so the time-integration fixes do not thermostate.
<LI>don't forget to time-integrate both cores and Drude particles.
<LI>pressure is time-integrated only once by using <I>nve</I> for Drude particles and <I>nph</I> for atoms/cores (or vice versa). Do not use <I>nph</I> for both.
<LI>contrary to the alternative thermostating using Nose-Hoover thermostat fix <I>npt</I> and fix <I>drude/transform</I>, the <I>fix_modify</I> command is not required here, because the fix <I>nph</I> computes the global pressure even if its group is <I>ATOMS</I>. This is what we want. If we thermostated <I>ATOMS</I> using <I>npt</I>, the pressure should be the global one, but the temperature should be only that of the cores. That's why the command <I>fix_modify</I> should be called in that case.
<LI>f_TEMP[1] and f_TEMP[2] contain the reduced temperatures of the cores/atoms and of the Drude particles (see below). They should be 300. and 1. on average here.
</UL>
<HR>
<P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
</P>
<P>No information about this fix is written to <A HREF = "restart.html">binary restart
files</A>. Because the state of the random number generator
is not saved in restart files, this means you cannot do "exact"
restarts with this fix, where the simulation continues on the same as
if no restart had taken place. However, in a statistical sense, a
restarted simulation should produce the same behavior.
</P>
<P>The <A HREF = "fix_modify.html">fix_modify</A> <I>temp</I> option is supported by this
fix. You can use it to assign a temperature <A HREF = "compute.html">compute</A>
you have defined to this fix which will be used in its thermostating
procedure, as described above. For consistency, the group used by the
compute should include the group of this fix and the Drude particles.
</P>
<P>This fix computes a global vector with 6 components which can be
accessed by various <A HREF = "Section_howto.html#howto_15">output commands</A>. The
meaning of the components are as follows:
</P>
<OL><LI>temperature of the centers of mass (temperature units)
<LI>temperature of the dipoles (temperature units)
<LI>number of degrees of freedom of the centers of mass
<LI>number of degrees of freedom of the dipoles
<LI>kinetic energy of the centers of mass (energy units)
<LI>kinetic energy of the dipoles (energy units)
</OL>
<P>This fix is not invoked during <A HREF = "minimize.html">energy minimization</A>.
</P>
<P><B>Restrictions:</B> none
</P>
<P><B>Related commands:</B>
</P>
<P><A HREF = "fix_langevin.html">fix langevin</A>,
<A HREF = "fix_drude_transform.html">fix drude/transform</A>,
<A HREF = "compute_temp_drude.html">compute temp/drude</A>,
<A HREF = "pair_thole.html">pair_style thole</A>
</P>
<P><B>Default:</B>
</P>
<P>The option defaults are zero = no.
</P>
<HR>
<A NAME = "Jiang"></A>
<P><B>(Jiang)</B> Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux, J
Phys Chem Lett, 2, 87-92 (2011).
</P>
</HTML>

253
doc/fix_langevin_drude.txt Normal file
View File

@ -0,0 +1,253 @@
<script type="text/javascript"
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
</script>
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix langevin/drude command :h3
[Syntax:]
fix ID group-ID langevin/drude Tcom damp_com seed_com Tdrude damp_drude seed_drude keyword values ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
langevin/drude = style name of this fix command :l
Tcom = desired temperature of the centers of mass (temperature units) :l
damp_com = damping parameter for the thermostat on centers of mass (time units) :l
seed_com = random number seed to use for white noise of the thermostat on centers of mass (positive integer) :l
Tdrude = desired temperature of the Drude oscillators (temperature units) :l
damp_drude = damping parameter for the thermostat on Drude oscillators (time units) :l
seed_drude = random number seed to use for white noise of the thermostat on Drude oscillators (positive integer) :l
zero or more keyword/value pairs may be appended :l
keyword = {zero} :l
{zero} value = {no} or {yes}
{no} = do not set total random force on centers of mass to zero
{yes} = set total random force on centers of mass to zero :pre
:ule
[Examples:]
fix 3 all langevin/drude 300.0 100.0 19377 1.0 20.0 83451
fix 1 all langevin/drude 298.15 100.0 19377 5.0 10.0 83451 zero yes :pre
[Description:]
Apply two Langevin thermostats as described in "(Jiang)"_#Jiang for
thermalizing the reduced degrees of freedom of Drude oscillators.
This link describes how to use the "thermalized Drude oscillator
model"_tutorial_drude.html in LAMMPS and polarizable models in LAMMPS
are discussed in "this Section"_Section_howto.html#howto_25.
Drude oscillators are a way to simulate polarizables atoms, by
splitting them into a core and a Drude particle bound by a harmonic
bond. The thermalization works by transforming the particles degrees
of freedom by these equations. In these equations upper case denotes
atomic or center of mass values and lower case denotes Drude particle
or dipole values. Primes denote the transformed (reduced) values,
while bare letters denote the original values.
Velocities:
\begin\{equation\} V' = \frac \{M\, V + m\, v\} \{M'\} \end\{equation\}
\begin\{equation\} v' = v - V \end\{equation\}
Masses:
\begin\{equation\} M' = M + m \end\{equation\}
\begin\{equation\} m' = \frac \{M\, m \} \{M'\} \end\{equation\}
The Langevin forces are computed as
\begin\{equation\} F' = - \frac \{M'\} \{\mathtt\{damp\_com\}\}\, V' + F_r' \end\{equation\}
\begin\{equation\} f' = - \frac \{m'\} \{\mathtt\{damp\_drude\}\}\, v' + f_r' \end\{equation\}
\(F_r'\) is a random force proportional to
\(\sqrt \{ \frac \{2\, k_B \mathtt\{Tcom\}\, m'\}
\{\mathrm dt\, \mathtt\{damp\_com\} \}
\} \). :b
\(f_r'\) is a random force proportional to
\(\sqrt \{ \frac \{2\, k_B \mathtt\{Tdrude\}\, m'\}
\{\mathrm dt\, \mathtt\{damp\_drude\} \}
\} \). :b
Then the real forces acting on the particles are computed from the inverse
transform:
\begin\{equation\} F = \frac M \{M'\}\, F' - f' \end\{equation\}
\begin\{equation\} f = \frac m \{M'\}\, F' + f' \end\{equation\}
For Drude pairs (core + electron), the center of mass and the dipole
are thermostated if (and only if) the core atom is in the specified
group.
Note that the thermostat effect of this fix is applied to only the
translational degrees of freedom of the particles, which is an
important consideration if finite-size particles, which have
rotational degrees of freedom, are being thermostated. The
translational degrees of freedom can also have a bias velocity removed
from them before thermostating takes place; see the description below.
IMPORTANT NOTE: Like the "fix langevin"_fix_langevin.html command,
this fix does NOT perform time integration. It only modifies forces to
effect thermostating. Thus you must use a separate time integration
fix, like "fix nve"_fix_nve.html or "fix nph"_fix_nh.html to actually
update the velocities and positions of atoms using the modified
forces. Likewise, this fix should not normally be used on atoms that
also have their temperature controlled by another fix - e.g. by "fix
nvt"_fix_nh.html or "fix temp/rescale"_fix_temp_rescale.html commands.
See "this howto section"_Section_howto.html#howto_16 of the manual for
a discussion of different ways to compute temperature and perform
thermostating.
:line
This fix requires each atom know whether it is a Drude particle or
not. You must therefore use the "fix drude"_fix_drude.html command to
specify the Drude status of each atom type.
IMPORTANT NOTE: only the Drude core atoms need to be in the group
specified for this fix. A Drude electron will be transformed together
with its cores even if it is not itself in the group. It is safe to
include Drude electrons or non-polarizable atoms in the group. The
non-polarizable atoms will simply not be thermostatted as if they had
a massless Drude partner (electron).
IMPORTANT NOTE: Ghost atoms need to know their velocity for this fix
to act correctly. You must use the "comm_modify"_comm_modify.html
command to enable this, e.g.
comm_modify vel yes :pre
:line
{Tcom} is the target temperature of the centers of mass, which would
be used to thermostate the non-polarizable atoms. {Tdrude} is the
(normally low) target temperature of the core-Drude particle pairs
(dipoles). {Tcom} and {Tdrude} can be specified as an equal-style
"variable"_variable.html. If the value is a variable, it should be
specified as v_name, where name is the variable name. In this case,
the variable will be evaluated each timestep, and its value used to
determine the target temperature.
Equal-style variables can specify formulas with various mathematical
functions, and include "thermo_style"_thermo_style.html command
keywords for the simulation box parameters and timestep and elapsed
time. Thus it is easy to specify a time-dependent temperature.
Like other fixes that perform thermostating, this fix can be used with
"compute commands"_compute.html that remove a "bias" from the atom
velocities. E.g. removing the center-of-mass velocity from a group of
atoms. This is not done by default, but only if the
"fix_modify"_fix_modify.html command is used to assign a temperature
compute to this fix that includes such a bias term. See the doc pages
for individual "compute commands"_compute.html to determine which ones
include a bias. In this case, the thermostat works in the following
manner: bias is removed from each atom, thermostating is performed on
the remaining thermal degrees of freedom, and the bias is added back
in. NOTE: this feature has not been tested.
Note: The temperature thermostating the core-Drude particle pairs
should be chosen low enough, so as to mimic as closely as possible the
self-consistent minimization. It must however be high enough, so that
the dipoles can follow the local electric field exerted by the
neighbouring atoms. The optimal value probably depends on the
temperature of the centers of mass and on the mass of the Drude
particles.
{damp_com} is the characteristic time for reaching thermal equilibrium
of the centers of mass. For example, a value of 100.0 means to relax
the temperature of the centers of mass in a timespan of (roughly) 100
time units (tau or fmsec or psec - see the "units"_units.html
command). {damp_drude} is the characteristic time for reaching
thermal equilibrium of the dipoles. It is typically a few timesteps.
The number {seed_com} and {seed_drude} are positive integers. They set
the seeds of the Marsaglia random number generators used for
generating the random forces on centers of mass and on the
dipoles. Each processor uses the input seed to generate its own unique
seed and its own stream of random numbers. Thus the dynamics of the
system will not be identical on two runs on different numbers of
processors.
The keyword {zero} can be used to eliminate drift due to the
thermostat on centers of mass. Because the random forces on different
centers of mass are independent, they do not sum exactly to zero. As
a result, this fix applies a small random force to the entire system,
and the momentum of the total center of mass of the system undergoes a
slow random walk. If the keyword {zero} is set to {yes}, the total
random force on the centers of mass is set exactly to zero by
subtracting off an equal part of it from each center of mass in the
group. As a result, the total center of mass of a system with zero
initial momentum will not drift over time.
:line
Example for rigid bodies in the NPT ensemble:
comm_modify vel yes
fix TEMP all langevin/drude 300. 100. 1256 1. 20. 13977 zero yes
fix NPH ATOMS rigid/nph/small molecule iso 1. 1. 500.
fix NVE DRUDES nve
thermo_style custom step cpu etotal ke pe ebond ecoul elong press vol temp c_TATOM c_TDRUDE f_TEMP\[1\] f_TEMP\[2\] :pre
Comments:
Drude particles should not be in the rigid group, otherwise the Drude oscillators will be frozen and the system will lose its polarizability.
{zero yes} avoids a drift of the center of mass of the system, but is a bit slower.
use two different random seeds to avoid unphysical correlations.
temperature is controlled by the fix {langevin/drude}, so the time-integration fixes do not thermostate.
don't forget to time-integrate both cores and Drude particles.
pressure is time-integrated only once by using {nve} for Drude particles and {nph} for atoms/cores (or vice versa). Do not use {nph} for both.
contrary to the alternative thermostating using Nose-Hoover thermostat fix {npt} and fix {drude/transform}, the {fix_modify} command is not required here, because the fix {nph} computes the global pressure even if its group is {ATOMS}. This is what we want. If we thermostated {ATOMS} using {npt}, the pressure should be the global one, but the temperature should be only that of the cores. That's why the command {fix_modify} should be called in that case.
f_TEMP\[1\] and f_TEMP\[2\] contain the reduced temperatures of the cores/atoms and of the Drude particles (see below). They should be 300. and 1. on average here. :ul
:line
[Restart, fix_modify, output, run start/stop, minimize info:]
No information about this fix is written to "binary restart
files"_restart.html. Because the state of the random number generator
is not saved in restart files, this means you cannot do "exact"
restarts with this fix, where the simulation continues on the same as
if no restart had taken place. However, in a statistical sense, a
restarted simulation should produce the same behavior.
The "fix_modify"_fix_modify.html {temp} option is supported by this
fix. You can use it to assign a temperature "compute"_compute.html
you have defined to this fix which will be used in its thermostating
procedure, as described above. For consistency, the group used by the
compute should include the group of this fix and the Drude particles.
This fix computes a global vector with 6 components which can be
accessed by various "output commands"_Section_howto.html#howto_15. The
meaning of the components are as follows:
temperature of the centers of mass (temperature units)
temperature of the dipoles (temperature units)
number of degrees of freedom of the centers of mass
number of degrees of freedom of the dipoles
kinetic energy of the centers of mass (energy units)
kinetic energy of the dipoles (energy units) :ol
This fix is not invoked during "energy minimization"_minimize.html.
[Restrictions:] none
[Related commands:]
"fix langevin"_fix_langevin.html,
"fix drude/transform"_fix_drude_transform.html,
"compute temp/drude"_compute_temp_drude.html,
"pair_style thole"_pair_thole.html
[Default:]
The option defaults are zero = no.
:line
:link(Jiang)
[(Jiang)] Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux, J
Phys Chem Lett, 2, 87-92 (2011).

4
doc/pair_bop.html
View File

@ -40,7 +40,9 @@ transferability to different phases can approach that of quantum
mechanical methods. This potential is similar to the original BOP mechanical methods. This potential is similar to the original BOP
developed by Pettifor (<A HREF = "#Pettifor_1">Pettifor_1</A>, developed by Pettifor (<A HREF = "#Pettifor_1">Pettifor_1</A>,
<A HREF = "#Pettifor_2">Pettifor_2</A>, <A HREF = "#Pettifor_3">Pettifor_3</A>) and later updated <A HREF = "#Pettifor_2">Pettifor_2</A>, <A HREF = "#Pettifor_3">Pettifor_3</A>) and later updated
by Murdick, Zhou, and Ward (<A HREF = "#Murdick">Murdick</A>, <A HREF = "#Ward">Ward</A>). by Murdick, Zhou, and Ward (<A HREF = "#Murdick">Murdick</A>, <A HREF = "#Ward">Ward</A>). As of
summer 2015, BOP potential files for these systems are provide with
LAMMPS: AlCu, CCu, CdTe, CuH, GaAs.
</P> </P>
<P>The BOP potential consists of three terms: <P>The BOP potential consists of three terms:
</P> </P>

4
doc/pair_bop.txt
View File

@ -34,7 +34,9 @@ transferability to different phases can approach that of quantum
mechanical methods. This potential is similar to the original BOP mechanical methods. This potential is similar to the original BOP
developed by Pettifor ("Pettifor_1"_#Pettifor_1, developed by Pettifor ("Pettifor_1"_#Pettifor_1,
"Pettifor_2"_#Pettifor_2, "Pettifor_3"_#Pettifor_3) and later updated "Pettifor_2"_#Pettifor_2, "Pettifor_3"_#Pettifor_3) and later updated
by Murdick, Zhou, and Ward ("Murdick"_#Murdick, "Ward"_#Ward). by Murdick, Zhou, and Ward ("Murdick"_#Murdick, "Ward"_#Ward). As of
summer 2015, BOP potential files for these systems are provide with
LAMMPS: AlCu, CCu, CdTe, CuH, GaAs.
The BOP potential consists of three terms: The BOP potential consists of three terms:

15
doc/pair_born.html
View File

@ -17,6 +17,8 @@
</H3> </H3>
<H3>pair_style born/coul/long command <H3>pair_style born/coul/long command
</H3> </H3>
<H3>pair_style born/coul/long/cs command
</H3>
<H3>pair_style born/coul/long/cuda command <H3>pair_style born/coul/long/cuda command
</H3> </H3>
<H3>pair_style born/coul/long/gpu command <H3>pair_style born/coul/long/gpu command
@ -37,12 +39,12 @@
</P> </P>
<PRE>pair_style style args <PRE>pair_style style args
</PRE> </PRE>
<UL><LI>style = <I>born</I> or <I>born/coul/long</I> or <I>born/coul/msm</I> or <I>born/coul/wolf</I> <UL><LI>style = <I>born</I> or <I>born/coul/long</I> or <I>born/coul/long/cs</I> or <I>born/coul/msm</I> or <I>born/coul/wolf</I>
<LI>args = list of arguments for a particular style <LI>args = list of arguments for a particular style
</UL> </UL>
<PRE> <I>born</I> args = cutoff <PRE> <I>born</I> args = cutoff
cutoff = global cutoff for non-Coulombic interactions (distance units) cutoff = global cutoff for non-Coulombic interactions (distance units)
<I>born/coul/long</I> args = cutoff (cutoff2) <I>born/coul/long</I> or <I>born/coul/long/cs</I> args = cutoff (cutoff2)
cutoff = global cutoff for non-Coulombic (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for non-Coulombic (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units)
<I>born/coul/msm</I> args = cutoff (cutoff2) <I>born/coul/msm</I> args = cutoff (cutoff2)
@ -60,7 +62,9 @@ pair_coeff * * 6.08 0.317 2.340 24.18 11.51
pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51
</PRE> </PRE>
<PRE>pair_style born/coul/long 10.0 <PRE>pair_style born/coul/long 10.0
pair_style born/coul/long/cs 10.0
pair_style born/coul/long 10.0 8.0 pair_style born/coul/long 10.0 8.0
pair_style born/coul/long/cs 10.0 8.0
pair_coeff * * 6.08 0.317 2.340 24.18 11.51 pair_coeff * * 6.08 0.317 2.340 24.18 11.51
pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51
</PRE> </PRE>
@ -102,7 +106,12 @@ term.
<P>The <I>born/coul/wolf</I> style adds a Coulombic term as described for the <P>The <I>born/coul/wolf</I> style adds a Coulombic term as described for the
Wolf potential in the <A HREF = "pair_coul.html">coul/wolf</A> pair style. Wolf potential in the <A HREF = "pair_coul.html">coul/wolf</A> pair style.
</P> </P>
<P>Note that these potentials are related to the <A HREF = "pair_born.html">Buckingham <P>Style <I>born/coul/long/cs</I> is identical to <I>born/coul/long</I> except that
a term is added for the <A HREF = "Section_howto.html#howto_25">core/shell model</A>
to allow charges on core and shell particles to be separated by r =
0.0.
</P>
<P>Note that these potentials are related to the <A HREF = "pair_buck.html">Buckingham
potential</A>. potential</A>.
</P> </P>
<P>The following coefficients must be defined for each pair of atoms <P>The following coefficients must be defined for each pair of atoms

14
doc/pair_born.txt
View File

@ -11,6 +11,7 @@ pair_style born command :h3
pair_style born/omp command :h3 pair_style born/omp command :h3
pair_style born/gpu command :h3 pair_style born/gpu command :h3
pair_style born/coul/long command :h3 pair_style born/coul/long command :h3
pair_style born/coul/long/cs command :h3
pair_style born/coul/long/cuda command :h3 pair_style born/coul/long/cuda command :h3
pair_style born/coul/long/gpu command :h3 pair_style born/coul/long/gpu command :h3
pair_style born/coul/long/omp command :h3 pair_style born/coul/long/omp command :h3
@ -24,11 +25,11 @@ pair_style born/coul/wolf/omp command :h3
pair_style style args :pre pair_style style args :pre
style = {born} or {born/coul/long} or {born/coul/msm} or {born/coul/wolf} style = {born} or {born/coul/long} or {born/coul/long/cs} or {born/coul/msm} or {born/coul/wolf}
args = list of arguments for a particular style :ul args = list of arguments for a particular style :ul
{born} args = cutoff {born} args = cutoff
cutoff = global cutoff for non-Coulombic interactions (distance units) cutoff = global cutoff for non-Coulombic interactions (distance units)
{born/coul/long} args = cutoff (cutoff2) {born/coul/long} or {born/coul/long/cs} args = cutoff (cutoff2)
cutoff = global cutoff for non-Coulombic (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for non-Coulombic (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units)
{born/coul/msm} args = cutoff (cutoff2) {born/coul/msm} args = cutoff (cutoff2)
@ -46,7 +47,9 @@ pair_coeff * * 6.08 0.317 2.340 24.18 11.51
pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre
pair_style born/coul/long 10.0 pair_style born/coul/long 10.0
pair_style born/coul/long/cs 10.0
pair_style born/coul/long 10.0 8.0 pair_style born/coul/long 10.0 8.0
pair_style born/coul/long/cs 10.0 8.0
pair_coeff * * 6.08 0.317 2.340 24.18 11.51 pair_coeff * * 6.08 0.317 2.340 24.18 11.51
pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre
@ -88,8 +91,13 @@ term.
The {born/coul/wolf} style adds a Coulombic term as described for the The {born/coul/wolf} style adds a Coulombic term as described for the
Wolf potential in the "coul/wolf"_pair_coul.html pair style. Wolf potential in the "coul/wolf"_pair_coul.html pair style.
Style {born/coul/long/cs} is identical to {born/coul/long} except that
a term is added for the "core/shell model"_Section_howto.html#howto_25
to allow charges on core and shell particles to be separated by r =
0.0.
Note that these potentials are related to the "Buckingham Note that these potentials are related to the "Buckingham
potential"_pair_born.html. potential"_pair_buck.html.
The following coefficients must be defined for each pair of atoms The following coefficients must be defined for each pair of atoms
types via the "pair_coeff"_pair_coeff.html command as in the examples types via the "pair_coeff"_pair_coeff.html command as in the examples

20
doc/pair_buck.html
View File

@ -31,6 +31,8 @@
</H3> </H3>
<H3>pair_style buck/coul/long command <H3>pair_style buck/coul/long command
</H3> </H3>
<H3>pair_style buck/coul/long/cs command
</H3>
<H3>pair_style buck/coul/long/cuda command <H3>pair_style buck/coul/long/cuda command
</H3> </H3>
<H3>pair_style buck/coul/long/gpu command <H3>pair_style buck/coul/long/gpu command
@ -47,7 +49,7 @@
</P> </P>
<PRE>pair_style style args <PRE>pair_style style args
</PRE> </PRE>
<UL><LI>style = <I>buck</I> or <I>buck/coul/cut</I> or <I>buck/coul/long</I> or <I>buck/coul/msm</I> <UL><LI>style = <I>buck</I> or <I>buck/coul/cut</I> or <I>buck/coul/long</I> or <I>buck/coul/long/cs</I> or <I>buck/coul/msm</I>
<LI>args = list of arguments for a particular style <LI>args = list of arguments for a particular style
</UL> </UL>
<PRE> <I>buck</I> args = cutoff <PRE> <I>buck</I> args = cutoff
@ -55,7 +57,7 @@
<I>buck/coul/cut</I> args = cutoff (cutoff2) <I>buck/coul/cut</I> args = cutoff (cutoff2)
cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units)
<I>buck/coul/long</I> args = cutoff (cutoff2) <I>buck/coul/long</I> or <I>buck/coul/long/cs</I> args = cutoff (cutoff2)
cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units)
<I>buck/coul/msm</I> args = cutoff (cutoff2) <I>buck/coul/msm</I> args = cutoff (cutoff2)
@ -75,7 +77,9 @@ pair_coeff 1 1 100.0 1.5 200.0 9.0
pair_coeff 1 1 100.0 1.5 200.0 9.0 8.0 pair_coeff 1 1 100.0 1.5 200.0 9.0 8.0
</PRE> </PRE>
<PRE>pair_style buck/coul/long 10.0 <PRE>pair_style buck/coul/long 10.0
pair_style buck/coul/long/cs 10.0
pair_style buck/coul/long 10.0 8.0 pair_style buck/coul/long 10.0 8.0
pair_style buck/coul/long/cs 10.0 8.0
pair_coeff * * 100.0 1.5 200.0 pair_coeff * * 100.0 1.5 200.0
pair_coeff 1 1 100.0 1.5 200.0 9.0 pair_coeff 1 1 100.0 1.5 200.0 9.0
</PRE> </PRE>
@ -109,6 +113,11 @@ A,C and Coulombic terms. If two cutoffs are specified, the first is
used as the cutoff for the A,C terms, and the second is the cutoff for used as the cutoff for the A,C terms, and the second is the cutoff for
the Coulombic term. the Coulombic term.
</P> </P>
<P>Style <I>buck/coul/long/cs</I> is identical to <I>buck/coul/long</I> except that
a term is added for the <A HREF = "Section_howto.html#howto_25">core/shell model</A>
to allow charges on core and shell particles to be separated by r =
0.0.
</P>
<P>Note that these potentials are related to the <A HREF = "pair_born.html">Born-Mayer-Huggins <P>Note that these potentials are related to the <A HREF = "pair_born.html">Born-Mayer-Huggins
potential</A>. potential</A>.
</P> </P>
@ -179,7 +188,7 @@ I,J pairs must be specified explicitly.
for the energy of the exp() and 1/r^6 portion of the pair interaction. for the energy of the exp() and 1/r^6 portion of the pair interaction.
</P> </P>
<P>The <I>buck/coul/long</I> pair style supports the <P>The <I>buck/coul/long</I> pair style supports the
<A HREF = "pair_modify.html">pair_modify</A> table option ti tabulate the <A HREF = "pair_modify.html">pair_modify</A> table option to tabulate the
short-range portion of the long-range Coulombic interaction. short-range portion of the long-range Coulombic interaction.
</P> </P>
<P>These styles support the pair_modify tail option for adding long-range <P>These styles support the pair_modify tail option for adding long-range
@ -196,8 +205,9 @@ respa</A> command. They do not support the <I>inner</I>,
</P> </P>
<P><B>Restrictions:</B> <P><B>Restrictions:</B>
</P> </P>
<P>The <I>buck/coul/long</I> style is part of the KSPACE package. It is only <P>The <I>buck/coul/long</I> style is part of the KSPACE package. The
enabled if LAMMPS was built with that package (which it is by <I>buck/coul/long/cs</I> style is part of the CORESHELL package. They are
only enabled if LAMMPS was built with that package (which it is by
default). See the <A HREF = "Section_start.html#start_3">Making LAMMPS</A> section default). See the <A HREF = "Section_start.html#start_3">Making LAMMPS</A> section
for more info. for more info.
</P> </P>

19
doc/pair_buck.txt
View File

@ -17,6 +17,7 @@ pair_style buck/coul/cut/gpu command :h3
pair_style buck/coul/cut/kk command :h3 pair_style buck/coul/cut/kk command :h3
pair_style buck/coul/cut/omp command :h3 pair_style buck/coul/cut/omp command :h3
pair_style buck/coul/long command :h3 pair_style buck/coul/long command :h3
pair_style buck/coul/long/cs command :h3
pair_style buck/coul/long/cuda command :h3 pair_style buck/coul/long/cuda command :h3
pair_style buck/coul/long/gpu command :h3 pair_style buck/coul/long/gpu command :h3
pair_style buck/coul/long/kk command :h3 pair_style buck/coul/long/kk command :h3
@ -28,14 +29,14 @@ pair_style buck/coul/msm/omp command :h3
pair_style style args :pre pair_style style args :pre
style = {buck} or {buck/coul/cut} or {buck/coul/long} or {buck/coul/msm} style = {buck} or {buck/coul/cut} or {buck/coul/long} or {buck/coul/long/cs} or {buck/coul/msm}
args = list of arguments for a particular style :ul args = list of arguments for a particular style :ul
{buck} args = cutoff {buck} args = cutoff
cutoff = global cutoff for Buckingham interactions (distance units) cutoff = global cutoff for Buckingham interactions (distance units)
{buck/coul/cut} args = cutoff (cutoff2) {buck/coul/cut} args = cutoff (cutoff2)
cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units)
{buck/coul/long} args = cutoff (cutoff2) {buck/coul/long} or {buck/coul/long/cs} args = cutoff (cutoff2)
cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units) cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units)
{buck/coul/msm} args = cutoff (cutoff2) {buck/coul/msm} args = cutoff (cutoff2)
@ -55,7 +56,9 @@ pair_coeff 1 1 100.0 1.5 200.0 9.0
pair_coeff 1 1 100.0 1.5 200.0 9.0 8.0 :pre pair_coeff 1 1 100.0 1.5 200.0 9.0 8.0 :pre
pair_style buck/coul/long 10.0 pair_style buck/coul/long 10.0
pair_style buck/coul/long/cs 10.0
pair_style buck/coul/long 10.0 8.0 pair_style buck/coul/long 10.0 8.0
pair_style buck/coul/long/cs 10.0 8.0
pair_coeff * * 100.0 1.5 200.0 pair_coeff * * 100.0 1.5 200.0
pair_coeff 1 1 100.0 1.5 200.0 9.0 :pre pair_coeff 1 1 100.0 1.5 200.0 9.0 :pre
@ -89,6 +92,11 @@ A,C and Coulombic terms. If two cutoffs are specified, the first is
used as the cutoff for the A,C terms, and the second is the cutoff for used as the cutoff for the A,C terms, and the second is the cutoff for
the Coulombic term. the Coulombic term.
Style {buck/coul/long/cs} is identical to {buck/coul/long} except that
a term is added for the "core/shell model"_Section_howto.html#howto_25
to allow charges on core and shell particles to be separated by r =
0.0.
Note that these potentials are related to the "Born-Mayer-Huggins Note that these potentials are related to the "Born-Mayer-Huggins
potential"_pair_born.html. potential"_pair_born.html.
@ -159,7 +167,7 @@ These styles support the "pair_modify"_pair_modify.html shift option
for the energy of the exp() and 1/r^6 portion of the pair interaction. for the energy of the exp() and 1/r^6 portion of the pair interaction.
The {buck/coul/long} pair style supports the The {buck/coul/long} pair style supports the
"pair_modify"_pair_modify.html table option ti tabulate the "pair_modify"_pair_modify.html table option to tabulate the
short-range portion of the long-range Coulombic interaction. short-range portion of the long-range Coulombic interaction.
These styles support the pair_modify tail option for adding long-range These styles support the pair_modify tail option for adding long-range
@ -176,8 +184,9 @@ respa"_run_style.html command. They do not support the {inner},
[Restrictions:] [Restrictions:]
The {buck/coul/long} style is part of the KSPACE package. It is only The {buck/coul/long} style is part of the KSPACE package. The
enabled if LAMMPS was built with that package (which it is by {buck/coul/long/cs} style is part of the CORESHELL package. They are
only enabled if LAMMPS was built with that package (which it is by
default). See the "Making LAMMPS"_Section_start.html#start_3 section default). See the "Making LAMMPS"_Section_start.html#start_3 section
for more info. for more info.

21
doc/pair_coul.html
View File

@ -35,6 +35,8 @@
</H3> </H3>
<H3>pair_style coul/long command <H3>pair_style coul/long command
</H3> </H3>
<H3>pair_style coul/long/cs command
</H3>
<H3>pair_style coul/long/omp command <H3>pair_style coul/long/omp command
</H3> </H3>
<H3>pair_style coul/long/gpu command <H3>pair_style coul/long/gpu command
@ -67,6 +69,7 @@
pair_style coul/debye kappa cutoff pair_style coul/debye kappa cutoff
pair_style coul/dsf alpha cutoff pair_style coul/dsf alpha cutoff
pair_style coul/long cutoff pair_style coul/long cutoff
pair_style coul/long/cs cutoff
pair_style coul/long/gpu cutoff pair_style coul/long/gpu cutoff
pair_style coul/wolf alpha cutoff pair_style coul/wolf alpha cutoff
pair_style coul/streitz cutoff keyword alpha pair_style coul/streitz cutoff keyword alpha
@ -91,6 +94,7 @@ pair_coeff 2 2 3.5
pair_coeff * * pair_coeff * *
</PRE> </PRE>
<PRE>pair_style coul/long 10.0 <PRE>pair_style coul/long 10.0
pair_style coul/long/cs 10.0
pair_coeff * * pair_coeff * *
</PRE> </PRE>
<PRE>pair_style coul/msm 10.0 <PRE>pair_style coul/msm 10.0
@ -225,6 +229,10 @@ option. The Coulombic cutoff specified for this style means that
pairwise interactions within this distance are computed directly; pairwise interactions within this distance are computed directly;
interactions outside that distance are computed in reciprocal space. interactions outside that distance are computed in reciprocal space.
</P> </P>
<P>Style <I>coul/long/cs</I> is identical to <I>coul/long</I> except that a term is
added for the <A HREF = "Section_howto.html#howto_25">core/shell model</A> to allow
charges on core and shell particles to be separated by r = 0.0.
</P>
<P>Styles <I>tip4p/cut</I> and <I>tip4p/long</I> implement the coulomb part of <P>Styles <I>tip4p/cut</I> and <I>tip4p/long</I> implement the coulomb part of
the TIP4P water model of <A HREF = "#Jorgensen">(Jorgensen)</A>, which introduces the TIP4P water model of <A HREF = "#Jorgensen">(Jorgensen)</A>, which introduces
a massless site located a short distance away from the oxygen atom a massless site located a short distance away from the oxygen atom
@ -333,16 +341,15 @@ to be specified in an input script that reads a restart file.
<P><B>Restrictions:</B> <P><B>Restrictions:</B>
</P> </P>
<P>The <I>coul/long</I>, <I>coul/msm</I> and <I>tip4p/long</I> styles are part of the <P>The <I>coul/long</I>, <I>coul/msm</I> and <I>tip4p/long</I> styles are part of the
KSPACE package. They are only enabled if LAMMPS was built with that KSPACE package. The <I>coul/long/cs</I> style is part of the CORESHELL
package (which it is by default). package. They are only enabled if LAMMPS was built with that package
See the <A HREF = "Section_start.html#start_3">Making LAMMPS</A> section (which it is by default). See the <A HREF = "Section_start.html#start_3">Making
for more info. LAMMPS</A> section for more info.
</P> </P>
<P><B>Related commands:</B> <P><B>Related commands:</B>
</P> </P>
<P><A HREF = "pair_coeff.html">pair_coeff</A>, <A HREF = "pair_hybrid.html">pair_style <P><A HREF = "pair_coeff.html">pair_coeff</A>, <A HREF = "pair_hybrid.html">pair_style,
hybrid/overlay</A> hybrid/overlay</A>, <A HREF = "kspace_style.html">kspace_style</A>
<A HREF = "kspace_style.html">kspace_style</A>
</P> </P>
<P><B>Default:</B> none <P><B>Default:</B> none
</P> </P>

20
doc/pair_coul.txt
View File

@ -19,6 +19,7 @@ pair_style coul/dsf/gpu command :h3
pair_style coul/dsf/kk command :h3 pair_style coul/dsf/kk command :h3
pair_style coul/dsf/omp command :h3 pair_style coul/dsf/omp command :h3
pair_style coul/long command :h3 pair_style coul/long command :h3
pair_style coul/long/cs command :h3
pair_style coul/long/omp command :h3 pair_style coul/long/omp command :h3
pair_style coul/long/gpu command :h3 pair_style coul/long/gpu command :h3
pair_style coul/long/kk command :h3 pair_style coul/long/kk command :h3
@ -39,6 +40,7 @@ pair_style coul/cut cutoff
pair_style coul/debye kappa cutoff pair_style coul/debye kappa cutoff
pair_style coul/dsf alpha cutoff pair_style coul/dsf alpha cutoff
pair_style coul/long cutoff pair_style coul/long cutoff
pair_style coul/long/cs cutoff
pair_style coul/long/gpu cutoff pair_style coul/long/gpu cutoff
pair_style coul/wolf alpha cutoff pair_style coul/wolf alpha cutoff
pair_style coul/streitz cutoff keyword alpha pair_style coul/streitz cutoff keyword alpha
@ -63,6 +65,7 @@ pair_style coul/dsf 0.05 10.0
pair_coeff * * :pre pair_coeff * * :pre
pair_style coul/long 10.0 pair_style coul/long 10.0
pair_style coul/long/cs 10.0
pair_coeff * * :pre pair_coeff * * :pre
pair_style coul/msm 10.0 pair_style coul/msm 10.0
@ -197,6 +200,10 @@ option. The Coulombic cutoff specified for this style means that
pairwise interactions within this distance are computed directly; pairwise interactions within this distance are computed directly;
interactions outside that distance are computed in reciprocal space. interactions outside that distance are computed in reciprocal space.
Style {coul/long/cs} is identical to {coul/long} except that a term is
added for the "core/shell model"_Section_howto.html#howto_25 to allow
charges on core and shell particles to be separated by r = 0.0.
Styles {tip4p/cut} and {tip4p/long} implement the coulomb part of Styles {tip4p/cut} and {tip4p/long} implement the coulomb part of
the TIP4P water model of "(Jorgensen)"_#Jorgensen, which introduces the TIP4P water model of "(Jorgensen)"_#Jorgensen, which introduces
a massless site located a short distance away from the oxygen atom a massless site located a short distance away from the oxygen atom
@ -305,16 +312,15 @@ This pair style can only be used via the {pair} keyword of the
[Restrictions:] [Restrictions:]
The {coul/long}, {coul/msm} and {tip4p/long} styles are part of the The {coul/long}, {coul/msm} and {tip4p/long} styles are part of the
KSPACE package. They are only enabled if LAMMPS was built with that KSPACE package. The {coul/long/cs} style is part of the CORESHELL
package (which it is by default). package. They are only enabled if LAMMPS was built with that package
See the "Making LAMMPS"_Section_start.html#start_3 section (which it is by default). See the "Making
for more info. LAMMPS"_Section_start.html#start_3 section for more info.
[Related commands:] [Related commands:]
"pair_coeff"_pair_coeff.html, "pair_style "pair_coeff"_pair_coeff.html, "pair_style,
hybrid/overlay"_pair_hybrid.html hybrid/overlay"_pair_hybrid.html, "kspace_style"_kspace_style.html
"kspace_style"_kspace_style.html
[Default:] none [Default:] none

464
doc/tutorial_drude.html Normal file
View File

@ -0,0 +1,464 @@
<HTML>
<script type="text/javascript"
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
</script>
<H3>Tutorial for Thermalized Drude oscillators in LAMMPS
</H3>
<P>This tutorial explains how to use <A HREF = "drude_oscillators.html">Drude
oscillators</A> in LAMMPS to simulate polarizable
systems. As an illustration, the input files for a simulation of 250
phenol molecules are documented. First of all, LAMMPS has to be
compiled with the USER-DRUDE package activated. Then, the data file
and input scripts have to be modified to include the Drude dipoles and
how to handle them.
</P>
<HR>
<P><B>Overview of Drude induced dipoles</B>
</P>
<P>Polarizable atoms acquire an induced electric dipole moment under the
action of an external electric field, for example the electric field
created by the surrounding particles. Drude oscillators represent
these dipoles by two fixed charges: the core (DC) and the Drude
particle (DP) bound by a harmonic potential. The Drude particle can be
thought of as the electron cloud whose center can be displaced from
the position of the the corresponding nucleus.
</P>
<P>The sum of the masses of a core-Drude pair should be the mass of the
initial (unsplit) atom, \(m_C + m_D = m\). The sum of their charges
should be the charge of the initial (unsplit) atom, \(q_C + q_D = q\).
A harmonic potential between the core and Drude partners should be
present, with force constant \(k_D\) and an equilibrium distance of
zero. The (half-)stiffness of the <A HREF = "bond_harmonic.html">harmonic bond</A>
\(K_D = k_D/2\) and the Drude charge \(q_D\) are related to the atom
polarizability \(\alpha\) by
</P>
<P>\begin{equation} K_D = \frac 1 2\, \frac {q_D^2} \alpha
\end{equation}
</P>
<P>Ideally, the mass of the Drude particle should be small, and the
stiffness of the harmonic bond should be large, so that the Drude
particle remains close ot the core. The values of Drude mass, Drude
charge, and force constant can be chosen following different
strategies, as in the following examples of polarizable force
fields.
</P>
<UL><LI><A HREF = "#Lamoureux">Lamoureux and Roux</A> suggest adopting a global
half-stiffness, \(K_D\) = 500 kcal/(mol &Aring;<sup>2</sup>) &mdash;
which corresponds to a force constant \(k_D\) = 4184 kJ/(mol
&Aring;<sup>2</sup>) &mdash; for all types of core-Drude bond, a
global mass \(m_D\) = 0.4 g/mol (or u) for all types of Drude
particle, and to calculate the Drude charges for individual atom types
from the atom polarizabilities using equation (1). This choice is
followed in the polarizable CHARMM force field.
<LI><A HREF = "#Schroeder">Schroeder and Steinhauser</A> suggest adopting a global
charge \(q_D\) = -1.0e and a global mass \(m_D\) = 0.1 g/mol (or u)
for all Drude particles, and to calculate the force constant for each
type of core-Drude bond from equation (1). The timesteps used by these
authors are between 0.5 and 2 fs, with the degrees of freedom of the
Drude oscillators kept cold at 1 K. In both these force fields
hydrogen atoms are treated as non-polarizable.
</UL>
<P>The motion of of the Drude particles can be calculated by minimizing
the energy of the induced dipoles at each timestep, by an interative,
self-consistent procedure. The Drude particles can be massless and
therefore do not contribute to the kinetic energy. However, the
relaxed method is computationall slow. An extended-lagrangian method
can be used to calculate the positions of the Drude particles, but
this requires them to have mass. It is important in this case to
decouple the degrees of freedom associated with the Drude oscillators
from those of the normal atoms. Thermalizing the Drude dipoles at
temperatures comparable to the rest of the simulation leads to several
problems (kinetic energy transfer, very short timestep, etc.), which
can be remediated by the "cold Drude" technique (<A HREF = "#Lamoureux">Lamoureux and
Roux</A>).
</P>
<P>Two closely related models are used to represent polarization through
"charges on a spring": the core-shell model and the Drude
model. Although the basic idea is the same, the core-shell model is
normally used for ionic/crystalline materials, whereas the Drude model
is normally used for molecular systems and fluid states. In ionic
crystals the symmetry around each ion and the distance between them
are such that the core-shell model is sufficiently stable. But to be
applicable to molecular/covalent systems the Drude model includes two
important features:
</P>
<OL><LI>The possibility to thermostat the additional degrees of freedom
associated with the induced dipoles at very low temperature, in terms
of the reduced coordinates of the Drude particles with respect to
their cores. This makes the trajectory close to that of relaxed
induced dipoles.
<LI>The Drude dipoles on covalently bonded atoms interact too strongly
due to the short distances, so an atom may capture the Drude particle
(shell) of a neighbor, or the induced dipoles within the same molecule
may align too much. To avoid this, damping at short of the
interactions between the point charges composing the induced dipole
can be done by <A HREF = "#Thole">Thole</A> functions.
</OL>
<HR>
<P><B>Preparation of the data file</B>
</P>
<P>The data file is similar to a standard LAMMPS data file for
<I>atom_style full</I>. The DPs and the <I>harmonic bonds</I> connecting them
to their DC should appear in the data file as normal atoms and bonds.
</P>
<P>You can use the <I>polarizer</I> tool (Python script distributed with the
USER-DRUDE package) to convert a non-polarizable data file (here
<I>data.102494.lmp</I>) to a polarizable data file (<I>data-p.lmp</I>)
</P>
<PRE>polarizer -q -f phenol.dff data.102494.lmp data-p.lmp
</PRE>
<P>This will automatically insert the new atoms and bonds.
The masses and charges of DCs and DPs are computed
from <I>phenol.dff</I>, as well as the DC-DP bond constants. The file
<I>phenol.dff</I> contains the polarizabilities of the atom types
and the mass of the Drude particles, for instance:
</P>
<PRE># units: kJ/mol, A, deg
# kforce is in the form k/2 r_D^2
# type m_D/u q_D/e k_D alpha/A3 thole
OH 0.4 -1.0 4184.0 0.63 0.67
CA 0.4 -1.0 4184.0 1.36 2.51
CAI 0.4 -1.0 4184.0 1.09 2.51
</PRE>
<P>The hydrogen atoms are absent from this file, so they will be treated
as non-polarizable atoms. In the non-polarizable data file
<I>data.102494.lmp</I>, atom names corresponding to the atom type numbers
have to be specified as comments at the end of lines of the <I>Masses</I>
section. You probably need to edit it to add these names. It should
look like
</P>
<PRE>Masses
</PRE>
<PRE>1 12.011 # CAI
2 12.011 # CA
3 15.999 # OH
4 1.008 # HA
5 1.008 # HO
</PRE>
<HR>
<P><B>Basic input file</B>
</P>
<P>The atom style should be set to (or derive from) <I>full</I>, so that you
can define atomic charges and molecular bonds, angles, dihedrals...
</P>
<P>The <I>polarizer</I> tool also outputs certain lines related to the input
script (the use of these lines will be explained below). In order for
LAMMPS to recognize that you are using Drude oscillators, you should
use the fix <I>drude</I>. The command is
</P>
<PRE>fix DRUDE all drude 1 1 1 0 0 2 2 2
</PRE>
<P>or, equivalently
</P>
<PRE>fix DRUDE all drude C C C N N D D D
</PRE>
<P>The 0, 1, 2 (or N, C, D) following the <I>drude</I> keyword have the
following meaning: There is one tag for each atom type. This tag is 1
(or C) for DCs, 2 (or D) for DPs and 0 (or N) for non-polarizable
atoms. Here the atom types 1 to 3 (C and O atoms) are DC, atom types
4 and 5 (H atoms) are non-polarizable and the atom types 6 to 8 are
the newly created DPs.
</P>
<P>By recognizing the fix <I>drude</I>, LAMMPS will find and store matching
DC-DP pairs and will treat DP as equivalent to their DC in the
<I>special bonds</I> relations. It may be necessary to extend the space
for storing such special relations. In this case extra space should
be reserved by using the <I>extra</I> keyword of the <I>special_bonds</I>
command. With our phenol, there is 1 more special neighbor for which
space is required. Otherwise LAMMPS crashes and gives the required
value.
</P>
<PRE>special_bonds lj/coul 0.0 0.0 0.5 extra 1
</PRE>
<P>Let us assume we want to run a simple NVT simulation at 300 K. Note
that Drude oscillators need to be thermalized at a low temperature in
order to approximate a self-consistent field (SCF), therefore it is not
possible to simulate an NVE ensemble with this package. Since dipoles
are approximated by a charged DC-DP pair, the <I>pair_style</I> must
include Coulomb interactions, for instance <I>lj/cut/coul/long</I> with
<I>kspace_style pppm</I>. For example, with a cutoff of 10. and a precision
1.e-4:
</P>
<PRE>pair_style lj/cut/coul/long 10.0
kspace_style pppm 1.0e-4
</PRE>
<P>As compared to the non-polarizable input file, <I>pair_coeff</I> lines need
to be added for the DPs. Since the DPs have no Lennard-Jones
interactions, their <I>epsilon</I> is 0. so the only <I>pair_coeff</I> line
that needs to be added is
</P>
<PRE>pair_coeff * 6* 0.0 0.0 # All-DPs
</PRE>
<P>Now for the thermalization, the simplest choice is to use the <A HREF = "fix_langevin_drude.html">fix
langevin/drude</A>.
</P>
<PRE>fix LANG all langevin/drude 300. 100 12435 1. 20 13977
</PRE>
<P>This applies a Langevin thermostat at temperature 300. to the centers
of mass of the DC-DP pairs, with relaxation time 100 and with random
seed 12345. This fix applies also a Langevin thermostat at temperature
1. to the relative motion of the DPs around their DCs, with relaxation
time 20 and random seed 13977. Only the DCs and non-polarizable
atoms need to be in this fix's group. LAMMPS will thermostate the DPs
together with their DC. For this, ghost atoms need to know their
velocities. Thus you need to add the following command:
</P>
<PRE>comm_modify vel yes
</PRE>
<P>In order to avoid that the center of mass of the whole system
drifts due to the random forces of the Langevin thermostat on DCs, you
can add the <I>zero yes</I> option at the end of the fix line.
</P>
<P>If the fix <I>shake</I> is used to constrain the C-H bonds, it should be
invoked after the fix <I>langevin/drude</I> for more accuracy.
</P>
<PRE>fix SHAKE ATOMS shake 0.0001 20 0 t 4 5
</PRE>
<P>IMPORTANT NOTE: The group of the fix <I>shake</I> must not include the DPs.
If the group <I>ATOMS</I> is defined by non-DPs atom types, you could use
</P>
<P>Since the fix <I>langevin/drude</I> does not perform time integration (just
modification of forces but no position/velocity updates), the fix
<I>nve</I> should be used in conjunction.
</P>
<PRE>fix NVE all nve
</PRE>
<P>Finally, do not forget to update the atom type elements if you use
them in a <I>dump_modify ... element ...</I> command, by adding the element
type of the DPs. Here for instance
</P>
<PRE>dump DUMP all custom 10 dump.lammpstrj id mol type element x y z ix iy iz
dump_modify DUMP element C C O H H D D D
</PRE>
<P>The input file should now be ready for use!
</P>
<P>You will notice that the global temperature <I>thermo_temp</I> computed by
LAMMPS is not 300. K as wanted. This is because LAMMPS treats DPs as
standard atoms in his default compute. If you want to output the
temperatures of the DC-DP pair centers of mass and of the DPs relative
to their DCs, you should use <I>thermo_style custom</I> with respectively
<I>f_LANG[1]</I> and <I>f_LANG[2]</I>. These should be close to 300. and
1. on average.
</P>
<PRE>thermo_style custom step temp f_LANG[1] f_LANG[2]
</PRE>
<HR>
<P><B>Thole screening</B>
</P>
<P>Dipolar interactions represented by point charges on springs may not
be stable, for example if the atomic polarizability is too high for
instance, a DP can escape from its DC and be captured by another DC,
which makes the force and energy diverge and the simulation
crash. Even without reaching this extreme case, the correlation
between nearby dipoles on the same molecule may be exagerated. Often,
special bond relations prevent bonded neighboring atoms to see the
charge of each other's DP, so that the problem does not always appear.
It is possible to use screened dipole dipole interactions by using the
<A HREF = "pair_thole.html"><I>pair_style thole</I></A>. This is implemented as a
correction to the Coulomb pair_styles, which dampens at short distance
the interactions between the charges representing the induced dipoles.
It is to be used as <I>hybrid/overlay</I> with any standard <I>coul</I> pair
style. In our example, we would use
</P>
<PRE>pair_style hybrid/overlay lj/cut/coul/long 10.0 thole 2.6 10.0
</PRE>
<P>This tells LAMMPS that we are using two pair_styles. The first one is
as above (<I>lj/cut/coul/long 10.0</I>). The second one is a <I>thole</I>
pair_style with default screening factor 2.6 (<A HREF = "#Noskov">Noskov</A>) and
cutoff 10.0.
</P>
<P>Since <I>hybrid/overlay</I> does not support mixing rules, the interaction
coefficients of all the pairs of atom types with i < j should be
explicitly defined. The output of the <I>polarizer</I> script can be used
to complete the <I>pair_coeff</I> section of the input file. In our
example, this will look like:
</P>
<PRE>pair_coeff 1 1 lj/cut/coul/long 0.0700 3.550
pair_coeff 1 2 lj/cut/coul/long 0.0700 3.550
pair_coeff 1 3 lj/cut/coul/long 0.1091 3.310
pair_coeff 1 4 lj/cut/coul/long 0.0458 2.985
pair_coeff 2 2 lj/cut/coul/long 0.0700 3.550
pair_coeff 2 3 lj/cut/coul/long 0.1091 3.310
pair_coeff 2 4 lj/cut/coul/long 0.0458 2.985
pair_coeff 3 3 lj/cut/coul/long 0.1700 3.070
pair_coeff 3 4 lj/cut/coul/long 0.0714 2.745
pair_coeff 4 4 lj/cut/coul/long 0.0300 2.420
pair_coeff * 5 lj/cut/coul/long 0.0000 0.000
pair_coeff * 6* lj/cut/coul/long 0.0000 0.000
pair_coeff 1 1 thole 1.090 2.510
pair_coeff 1 2 thole 1.218 2.510
pair_coeff 1 3 thole 0.829 1.590
pair_coeff 1 6 thole 1.090 2.510
pair_coeff 1 7 thole 1.218 2.510
pair_coeff 1 8 thole 0.829 1.590
pair_coeff 2 2 thole 1.360 2.510
pair_coeff 2 3 thole 0.926 1.590
pair_coeff 2 6 thole 1.218 2.510
pair_coeff 2 7 thole 1.360 2.510
pair_coeff 2 8 thole 0.926 1.590
pair_coeff 3 3 thole 0.630 0.670
pair_coeff 3 6 thole 0.829 1.590
pair_coeff 3 7 thole 0.926 1.590
pair_coeff 3 8 thole 0.630 0.670
pair_coeff 6 6 thole 1.090 2.510
pair_coeff 6 7 thole 1.218 2.510
pair_coeff 6 8 thole 0.829 1.590
pair_coeff 7 7 thole 1.360 2.510
pair_coeff 7 8 thole 0.926 1.590
pair_coeff 8 8 thole 0.630 0.670
</PRE>
<P>For the <I>thole</I> pair style the coefficients are
</P>
<OL><LI>the atom polarizability in units of cubic length
<LI>the screening factor of the Thole function (optional, default value specified by the pair_style command)
<LI>the cutoff (optional, default value defined by the pair_style command)
</OL>
<P>The special neighbors have charge-charge and charge-dipole
interactions screened by the <I>coul</I> factors of the <I>special_bonds</I>
command (0., 0., and 0.5 in the example above). Without using the
pair_style <I>thole</I>, dipole-dipole interactions are screened by the
same factor. By using the pair_style <I>thole</I>, dipole-dipole
interactions are screened by Thole's function, whatever their special
relationship (except within each DC-DP pair of course). Consider for
example 1-2 neighbors: using the pair_style <I>thole</I>, their dipoles
will see each other (despite the <I>coul</I> factor being 0.) and the
interactions between these dipoles will be damped by Thole's
function.
</P>
<HR>
<P><B>Thermostats and barostats</B>
</P>
<P>Using a Nose-Hoover barostat with the <I>langevin/drude</I> thermostat is
straightforward using fix <I>nph</I> instead of <I>nve</I>. For example:
</P>
<PRE>fix NPH all nph iso 1. 1. 500
</PRE>
<P>It is also possible to use a Nose-Hoover instead of a Langevin
thermostat. This requires to use <A HREF = "fix_drude_transform.html"><I>fix
drude/transform</I></A> just before and after the
time intergation fixes. The <I>fix drude/transform/direct</I> converts the
atomic masses, positions, velocities and forces into a reduced
representation, where the DCs transform into the centers of mass of
the DC-DP pairs and the DPs transform into their relative position
with respect to their DC. The <I>fix drude/transform/inverse</I> performs
the reverse transformation. For a NVT simulation, with the DCs and
atoms at 300 K and the DPs at 1 K relative to their DC one would use
</P>
<PRE>fix DIRECT all drude/transform/direct
fix NVT1 ATOMS nvt temp 300. 300. 100
fix NVT2 DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse
</PRE>
<P>For our phenol example, the groups would be defined as
</P>
<PRE>group ATOMS type 1 2 3 4 5 # DCs and non-polarizable atoms
group CORES type 1 2 3 # DCs
group DRUDES type 6 7 8 # DPs
</PRE>
<P>Note that with the fixes <I>drude/transform</I>, it is not required to
specify <I>comm_modify vel yes</I> because the fixes do it anyway (several
times and for the forces also). To avoid the flying ice cube artifact
<A HREF = "#Lamoureux">(Lamoureux)</A>, where the atoms progressively freeze and the
center of mass of the whole system drifts faster and faster, the <I>fix
momentum</I> can be used. For instance:
</P>
<PRE>fix MOMENTUM all momentum 100 linear 1 1 1
</PRE>
<P>It is a bit more tricky to run a NPT simulation with Nose-Hoover
barostat and thermostat. First, the volume should be integrated only
once. So the fix for DCs and atoms should be <I>npt</I> while the fix for
DPs should be <I>nvt</I> (or vice versa). Second, the <I>fix npt</I> computes a
global pressure and thus a global temperature whatever the fix group.
We do want the pressure to correspond to the whole system, but we want
the temperature to correspond to the fix group only. We must then use
the <I>fix_modify</I> command for this. In the end, the block of
instructions for thermostating and barostating will look like
</P>
<PRE>compute TATOMS ATOMS temp
fix DIRECT all drude/transform/direct
fix NPT ATOMS npt temp 300. 300. 100 iso 1. 1. 500
fix_modify NPT temp TATOMS press thermo_press
fix NVT DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse
</PRE>
<HR>
<P><B>Rigid bodies</B>
</P>
<P>You may want to simulate molecules as rigid bodies (but polarizable).
Common cases are water models such as <A HREF = "#SWM4-NDP">SWM4-NDP</A>, which is a
kind of polarizable TIP4P water. The rigid bodies and the DPs should
be integrated separately, even with the Langevin thermostat. Let us
review the different thermostats and ensemble combinations.
</P>
<P>NVT ensemble using Langevin thermostat:
</P>
<PRE>comm_modify vel yes
fix LANG all langevin/drude 300. 100 12435 1. 20 13977
fix RIGID ATOMS rigid/nve/small molecule
fix NVE DRUDES nve
</PRE>
<P>NVT ensemble using Nose-Hoover thermostat:
</P>
<PRE>fix DIRECT all drude/transform/direct
fix RIGID ATOMS rigid/nvt/small molecule temp 300. 300. 100
fix NVT DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse
</PRE>
<P>NPT ensemble with Langevin thermostat:
</P>
<PRE>comm_modify vel yes
fix LANG all langevin/drude 300. 100 12435 1. 20 13977
fix RIGID ATOMS rigid/nph/small molecule iso 1. 1. 500
fix NVE DRUDES nve
</PRE>
<P>NPT ensemble using Nose-Hoover thermostat:
</P>
<PRE>compute TATOM ATOMS temp
fix DIRECT all drude/transform/direct
fix RIGID ATOMS rigid/npt/small molecule temp 300. 300. 100 iso 1. 1. 500
fix_modify RIGID temp TATOM press thermo_press
fix NVT DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse
</PRE>
<HR>
<A NAME = "Lamoureux"></A>
<P><B>(Lamoureux)</B> Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003)
</P>
<A NAME = "Schroeder"></A>
<P><B>(Schroeder)</B> Schr&ouml;der and Steinhauser, J Chem Phys, 133,
154511 (2010).
</P>
<A NAME = "Jiang"></A>
<P><B>(Jiang)</B> Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux,
J Phys Chem Lett, 2, 87-92 (2011).
</P>
<A NAME = "Thole"></A>
<P><B>(Thole)</B> Chem Phys, 59, 341 (1981).
</P>
<A NAME = "Noskov"></A>
<P><B>(Noskov)</B> Noskov, Lamoureux and Roux, J Phys Chem B, 109, 6705 (2005).
</P>
<A NAME = "SWM4-NDP"></A>
<P><B>(SWM4-NDP)</B> Lamoureux, Harder, Vorobyov, Roux, MacKerell, Chem Phys
Let, 418, 245-249 (2006)
</P>
</HTML>

456
doc/tutorial_drude.txt Normal file
View File

@ -0,0 +1,456 @@
<script type="text/javascript"
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
</script>
Tutorial for Thermalized Drude oscillators in LAMMPS :h3
This tutorial explains how to use "Drude
oscillators"_drude_oscillators.html in LAMMPS to simulate polarizable
systems. As an illustration, the input files for a simulation of 250
phenol molecules are documented. First of all, LAMMPS has to be
compiled with the USER-DRUDE package activated. Then, the data file
and input scripts have to be modified to include the Drude dipoles and
how to handle them.
:line
[Overview of Drude induced dipoles]
Polarizable atoms acquire an induced electric dipole moment under the
action of an external electric field, for example the electric field
created by the surrounding particles. Drude oscillators represent
these dipoles by two fixed charges: the core (DC) and the Drude
particle (DP) bound by a harmonic potential. The Drude particle can be
thought of as the electron cloud whose center can be displaced from
the position of the the corresponding nucleus.
The sum of the masses of a core-Drude pair should be the mass of the
initial (unsplit) atom, \(m_C + m_D = m\). The sum of their charges
should be the charge of the initial (unsplit) atom, \(q_C + q_D = q\).
A harmonic potential between the core and Drude partners should be
present, with force constant \(k_D\) and an equilibrium distance of
zero. The (half-)stiffness of the "harmonic bond"_bond_harmonic.html
\(K_D = k_D/2\) and the Drude charge \(q_D\) are related to the atom
polarizability \(\alpha\) by
\begin\{equation\} K_D = \frac 1 2\, \frac \{q_D^2\} \alpha
\end\{equation\}
Ideally, the mass of the Drude particle should be small, and the
stiffness of the harmonic bond should be large, so that the Drude
particle remains close ot the core. The values of Drude mass, Drude
charge, and force constant can be chosen following different
strategies, as in the following examples of polarizable force
fields.
"Lamoureux and Roux"_#Lamoureux suggest adopting a global
half-stiffness, \(K_D\) = 500 kcal/(mol &Aring;<sup>2</sup>) &mdash;
which corresponds to a force constant \(k_D\) = 4184 kJ/(mol
&Aring;<sup>2</sup>) &mdash; for all types of core-Drude bond, a
global mass \(m_D\) = 0.4 g/mol (or u) for all types of Drude
particle, and to calculate the Drude charges for individual atom types
from the atom polarizabilities using equation (1). This choice is
followed in the polarizable CHARMM force field. :ulb,l
"Schroeder and Steinhauser"_#Schroeder suggest adopting a global
charge \(q_D\) = -1.0e and a global mass \(m_D\) = 0.1 g/mol (or u)
for all Drude particles, and to calculate the force constant for each
type of core-Drude bond from equation (1). The timesteps used by these
authors are between 0.5 and 2 fs, with the degrees of freedom of the
Drude oscillators kept cold at 1 K. In both these force fields
hydrogen atoms are treated as non-polarizable. :ule,l
The motion of of the Drude particles can be calculated by minimizing
the energy of the induced dipoles at each timestep, by an interative,
self-consistent procedure. The Drude particles can be massless and
therefore do not contribute to the kinetic energy. However, the
relaxed method is computationall slow. An extended-lagrangian method
can be used to calculate the positions of the Drude particles, but
this requires them to have mass. It is important in this case to
decouple the degrees of freedom associated with the Drude oscillators
from those of the normal atoms. Thermalizing the Drude dipoles at
temperatures comparable to the rest of the simulation leads to several
problems (kinetic energy transfer, very short timestep, etc.), which
can be remediated by the "cold Drude" technique ("Lamoureux and
Roux"_#Lamoureux).
Two closely related models are used to represent polarization through
"charges on a spring": the core-shell model and the Drude
model. Although the basic idea is the same, the core-shell model is
normally used for ionic/crystalline materials, whereas the Drude model
is normally used for molecular systems and fluid states. In ionic
crystals the symmetry around each ion and the distance between them
are such that the core-shell model is sufficiently stable. But to be
applicable to molecular/covalent systems the Drude model includes two
important features:
The possibility to thermostat the additional degrees of freedom
associated with the induced dipoles at very low temperature, in terms
of the reduced coordinates of the Drude particles with respect to
their cores. This makes the trajectory close to that of relaxed
induced dipoles. :olb,l
The Drude dipoles on covalently bonded atoms interact too strongly
due to the short distances, so an atom may capture the Drude particle
(shell) of a neighbor, or the induced dipoles within the same molecule
may align too much. To avoid this, damping at short of the
interactions between the point charges composing the induced dipole
can be done by "Thole"_#Thole functions. :ole,l
:line
[Preparation of the data file]
The data file is similar to a standard LAMMPS data file for
{atom_style full}. The DPs and the {harmonic bonds} connecting them
to their DC should appear in the data file as normal atoms and bonds.
You can use the {polarizer} tool (Python script distributed with the
USER-DRUDE package) to convert a non-polarizable data file (here
{data.102494.lmp}) to a polarizable data file ({data-p.lmp})
polarizer -q -f phenol.dff data.102494.lmp data-p.lmp :pre
This will automatically insert the new atoms and bonds.
The masses and charges of DCs and DPs are computed
from {phenol.dff}, as well as the DC-DP bond constants. The file
{phenol.dff} contains the polarizabilities of the atom types
and the mass of the Drude particles, for instance:
# units: kJ/mol, A, deg
# kforce is in the form k/2 r_D^2
# type m_D/u q_D/e k_D alpha/A3 thole
OH 0.4 -1.0 4184.0 0.63 0.67
CA 0.4 -1.0 4184.0 1.36 2.51
CAI 0.4 -1.0 4184.0 1.09 2.51 :pre
The hydrogen atoms are absent from this file, so they will be treated
as non-polarizable atoms. In the non-polarizable data file
{data.102494.lmp}, atom names corresponding to the atom type numbers
have to be specified as comments at the end of lines of the {Masses}
section. You probably need to edit it to add these names. It should
look like
Masses :pre
1 12.011 # CAI
2 12.011 # CA
3 15.999 # OH
4 1.008 # HA
5 1.008 # HO :pre
:line
[Basic input file]
The atom style should be set to (or derive from) {full}, so that you
can define atomic charges and molecular bonds, angles, dihedrals...
The {polarizer} tool also outputs certain lines related to the input
script (the use of these lines will be explained below). In order for
LAMMPS to recognize that you are using Drude oscillators, you should
use the fix {drude}. The command is
fix DRUDE all drude 1 1 1 0 0 2 2 2 :pre
or, equivalently
fix DRUDE all drude C C C N N D D D :pre
The 0, 1, 2 (or N, C, D) following the {drude} keyword have the
following meaning: There is one tag for each atom type. This tag is 1
(or C) for DCs, 2 (or D) for DPs and 0 (or N) for non-polarizable
atoms. Here the atom types 1 to 3 (C and O atoms) are DC, atom types
4 and 5 (H atoms) are non-polarizable and the atom types 6 to 8 are
the newly created DPs.
By recognizing the fix {drude}, LAMMPS will find and store matching
DC-DP pairs and will treat DP as equivalent to their DC in the
{special bonds} relations. It may be necessary to extend the space
for storing such special relations. In this case extra space should
be reserved by using the {extra} keyword of the {special_bonds}
command. With our phenol, there is 1 more special neighbor for which
space is required. Otherwise LAMMPS crashes and gives the required
value.
special_bonds lj/coul 0.0 0.0 0.5 extra 1 :pre
Let us assume we want to run a simple NVT simulation at 300 K. Note
that Drude oscillators need to be thermalized at a low temperature in
order to approximate a self-consistent field (SCF), therefore it is not
possible to simulate an NVE ensemble with this package. Since dipoles
are approximated by a charged DC-DP pair, the {pair_style} must
include Coulomb interactions, for instance {lj/cut/coul/long} with
{kspace_style pppm}. For example, with a cutoff of 10. and a precision
1.e-4:
pair_style lj/cut/coul/long 10.0
kspace_style pppm 1.0e-4 :pre
As compared to the non-polarizable input file, {pair_coeff} lines need
to be added for the DPs. Since the DPs have no Lennard-Jones
interactions, their {epsilon} is 0. so the only {pair_coeff} line
that needs to be added is
pair_coeff * 6* 0.0 0.0 # All-DPs :pre
Now for the thermalization, the simplest choice is to use the "fix
langevin/drude"_fix_langevin_drude.html.
fix LANG all langevin/drude 300. 100 12435 1. 20 13977 :pre
This applies a Langevin thermostat at temperature 300. to the centers
of mass of the DC-DP pairs, with relaxation time 100 and with random
seed 12345. This fix applies also a Langevin thermostat at temperature
1. to the relative motion of the DPs around their DCs, with relaxation
time 20 and random seed 13977. Only the DCs and non-polarizable
atoms need to be in this fix's group. LAMMPS will thermostate the DPs
together with their DC. For this, ghost atoms need to know their
velocities. Thus you need to add the following command:
comm_modify vel yes :pre
In order to avoid that the center of mass of the whole system
drifts due to the random forces of the Langevin thermostat on DCs, you
can add the {zero yes} option at the end of the fix line.
If the fix {shake} is used to constrain the C-H bonds, it should be
invoked after the fix {langevin/drude} for more accuracy.
fix SHAKE ATOMS shake 0.0001 20 0 t 4 5 :pre
IMPORTANT NOTE: The group of the fix {shake} must not include the DPs.
If the group {ATOMS} is defined by non-DPs atom types, you could use
Since the fix {langevin/drude} does not perform time integration (just
modification of forces but no position/velocity updates), the fix
{nve} should be used in conjunction.
fix NVE all nve :pre
Finally, do not forget to update the atom type elements if you use
them in a {dump_modify ... element ...} command, by adding the element
type of the DPs. Here for instance
dump DUMP all custom 10 dump.lammpstrj id mol type element x y z ix iy iz
dump_modify DUMP element C C O H H D D D :pre
The input file should now be ready for use!
You will notice that the global temperature {thermo_temp} computed by
LAMMPS is not 300. K as wanted. This is because LAMMPS treats DPs as
standard atoms in his default compute. If you want to output the
temperatures of the DC-DP pair centers of mass and of the DPs relative
to their DCs, you should use {thermo_style custom} with respectively
{f_LANG\[1\]} and {f_LANG\[2\]}. These should be close to 300. and
1. on average.
thermo_style custom step temp f_LANG\[1\] f_LANG\[2\] :pre
:line
[Thole screening]
Dipolar interactions represented by point charges on springs may not
be stable, for example if the atomic polarizability is too high for
instance, a DP can escape from its DC and be captured by another DC,
which makes the force and energy diverge and the simulation
crash. Even without reaching this extreme case, the correlation
between nearby dipoles on the same molecule may be exagerated. Often,
special bond relations prevent bonded neighboring atoms to see the
charge of each other's DP, so that the problem does not always appear.
It is possible to use screened dipole dipole interactions by using the
"{pair_style thole}"_pair_thole.html. This is implemented as a
correction to the Coulomb pair_styles, which dampens at short distance
the interactions between the charges representing the induced dipoles.
It is to be used as {hybrid/overlay} with any standard {coul} pair
style. In our example, we would use
pair_style hybrid/overlay lj/cut/coul/long 10.0 thole 2.6 10.0 :pre
This tells LAMMPS that we are using two pair_styles. The first one is
as above ({lj/cut/coul/long 10.0}). The second one is a {thole}
pair_style with default screening factor 2.6 ("Noskov"_#Noskov) and
cutoff 10.0.
Since {hybrid/overlay} does not support mixing rules, the interaction
coefficients of all the pairs of atom types with i < j should be
explicitly defined. The output of the {polarizer} script can be used
to complete the {pair_coeff} section of the input file. In our
example, this will look like:
pair_coeff 1 1 lj/cut/coul/long 0.0700 3.550
pair_coeff 1 2 lj/cut/coul/long 0.0700 3.550
pair_coeff 1 3 lj/cut/coul/long 0.1091 3.310
pair_coeff 1 4 lj/cut/coul/long 0.0458 2.985
pair_coeff 2 2 lj/cut/coul/long 0.0700 3.550
pair_coeff 2 3 lj/cut/coul/long 0.1091 3.310
pair_coeff 2 4 lj/cut/coul/long 0.0458 2.985
pair_coeff 3 3 lj/cut/coul/long 0.1700 3.070
pair_coeff 3 4 lj/cut/coul/long 0.0714 2.745
pair_coeff 4 4 lj/cut/coul/long 0.0300 2.420
pair_coeff * 5 lj/cut/coul/long 0.0000 0.000
pair_coeff * 6* lj/cut/coul/long 0.0000 0.000
pair_coeff 1 1 thole 1.090 2.510
pair_coeff 1 2 thole 1.218 2.510
pair_coeff 1 3 thole 0.829 1.590
pair_coeff 1 6 thole 1.090 2.510
pair_coeff 1 7 thole 1.218 2.510
pair_coeff 1 8 thole 0.829 1.590
pair_coeff 2 2 thole 1.360 2.510
pair_coeff 2 3 thole 0.926 1.590
pair_coeff 2 6 thole 1.218 2.510
pair_coeff 2 7 thole 1.360 2.510
pair_coeff 2 8 thole 0.926 1.590
pair_coeff 3 3 thole 0.630 0.670
pair_coeff 3 6 thole 0.829 1.590
pair_coeff 3 7 thole 0.926 1.590
pair_coeff 3 8 thole 0.630 0.670
pair_coeff 6 6 thole 1.090 2.510
pair_coeff 6 7 thole 1.218 2.510
pair_coeff 6 8 thole 0.829 1.590
pair_coeff 7 7 thole 1.360 2.510
pair_coeff 7 8 thole 0.926 1.590
pair_coeff 8 8 thole 0.630 0.670 :pre
For the {thole} pair style the coefficients are
the atom polarizability in units of cubic length
the screening factor of the Thole function (optional, default value specified by the pair_style command)
the cutoff (optional, default value defined by the pair_style command) :ol
The special neighbors have charge-charge and charge-dipole
interactions screened by the {coul} factors of the {special_bonds}
command (0., 0., and 0.5 in the example above). Without using the
pair_style {thole}, dipole-dipole interactions are screened by the
same factor. By using the pair_style {thole}, dipole-dipole
interactions are screened by Thole's function, whatever their special
relationship (except within each DC-DP pair of course). Consider for
example 1-2 neighbors: using the pair_style {thole}, their dipoles
will see each other (despite the {coul} factor being 0.) and the
interactions between these dipoles will be damped by Thole's
function.
:line
[Thermostats and barostats]
Using a Nose-Hoover barostat with the {langevin/drude} thermostat is
straightforward using fix {nph} instead of {nve}. For example:
fix NPH all nph iso 1. 1. 500 :pre
It is also possible to use a Nose-Hoover instead of a Langevin
thermostat. This requires to use "{fix
drude/transform}"_fix_drude_transform.html just before and after the
time intergation fixes. The {fix drude/transform/direct} converts the
atomic masses, positions, velocities and forces into a reduced
representation, where the DCs transform into the centers of mass of
the DC-DP pairs and the DPs transform into their relative position
with respect to their DC. The {fix drude/transform/inverse} performs
the reverse transformation. For a NVT simulation, with the DCs and
atoms at 300 K and the DPs at 1 K relative to their DC one would use
fix DIRECT all drude/transform/direct
fix NVT1 ATOMS nvt temp 300. 300. 100
fix NVT2 DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse :pre
For our phenol example, the groups would be defined as
group ATOMS type 1 2 3 4 5 # DCs and non-polarizable atoms
group CORES type 1 2 3 # DCs
group DRUDES type 6 7 8 # DPs :pre
Note that with the fixes {drude/transform}, it is not required to
specify {comm_modify vel yes} because the fixes do it anyway (several
times and for the forces also). To avoid the flying ice cube artifact
"(Lamoureux)"_#Lamoureux, where the atoms progressively freeze and the
center of mass of the whole system drifts faster and faster, the {fix
momentum} can be used. For instance:
fix MOMENTUM all momentum 100 linear 1 1 1 :pre
It is a bit more tricky to run a NPT simulation with Nose-Hoover
barostat and thermostat. First, the volume should be integrated only
once. So the fix for DCs and atoms should be {npt} while the fix for
DPs should be {nvt} (or vice versa). Second, the {fix npt} computes a
global pressure and thus a global temperature whatever the fix group.
We do want the pressure to correspond to the whole system, but we want
the temperature to correspond to the fix group only. We must then use
the {fix_modify} command for this. In the end, the block of
instructions for thermostating and barostating will look like
compute TATOMS ATOMS temp
fix DIRECT all drude/transform/direct
fix NPT ATOMS npt temp 300. 300. 100 iso 1. 1. 500
fix_modify NPT temp TATOMS press thermo_press
fix NVT DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse :pre
:line
[Rigid bodies]
You may want to simulate molecules as rigid bodies (but polarizable).
Common cases are water models such as "SWM4-NDP"_#SWM4-NDP, which is a
kind of polarizable TIP4P water. The rigid bodies and the DPs should
be integrated separately, even with the Langevin thermostat. Let us
review the different thermostats and ensemble combinations.
NVT ensemble using Langevin thermostat:
comm_modify vel yes
fix LANG all langevin/drude 300. 100 12435 1. 20 13977
fix RIGID ATOMS rigid/nve/small molecule
fix NVE DRUDES nve :pre
NVT ensemble using Nose-Hoover thermostat:
fix DIRECT all drude/transform/direct
fix RIGID ATOMS rigid/nvt/small molecule temp 300. 300. 100
fix NVT DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse :pre
NPT ensemble with Langevin thermostat:
comm_modify vel yes
fix LANG all langevin/drude 300. 100 12435 1. 20 13977
fix RIGID ATOMS rigid/nph/small molecule iso 1. 1. 500
fix NVE DRUDES nve :pre
NPT ensemble using Nose-Hoover thermostat:
compute TATOM ATOMS temp
fix DIRECT all drude/transform/direct
fix RIGID ATOMS rigid/npt/small molecule temp 300. 300. 100 iso 1. 1. 500
fix_modify RIGID temp TATOM press thermo_press
fix NVT DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse :pre
:line
:link(Lamoureux)
[(Lamoureux)] Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003)
:link(Schroeder)
[(Schroeder)] Schr&ouml;der and Steinhauser, J Chem Phys, 133,
154511 (2010).
:link(Jiang)
[(Jiang)] Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux,
J Phys Chem Lett, 2, 87-92 (2011).
:link(Thole)
[(Thole)] Chem Phys, 59, 341 (1981).
:link(Noskov)
[(Noskov)] Noskov, Lamoureux and Roux, J Phys Chem B, 109, 6705 (2005).
:link(SWM4-NDP)
[(SWM4-NDP)] Lamoureux, Harder, Vorobyov, Roux, MacKerell, Chem Phys
Let, 418, 245-249 (2006)