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@13464 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2015-05-15 16:33:17 +00:00
parent 2284dcc452
commit c3fd41010d
23 changed files with 10777 additions and 20440 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
examples/USER/fep/CH4hyd/fdti01/fdti01.lmp
View File

@ -1,22 +1,22 @@
# Time-averaged data for fix fFEP # Time-averaged data for fix FEP
# TimeStep c_cFEP[1] c_cFEP[2] # TimeStep c_FEP[1] c_FEP[2]
125000 0.000490778 0.999186 100000 0.000289416 0.999521
250000 0.00619786 0.989679 200000 0.00590502 0.990163
375000 0.0128045 0.978813 300000 0.0115179 0.980934
500000 0.0203102 0.96667 400000 0.0216052 0.96457
625000 0.0272161 0.955808 500000 0.0222451 0.963797
750000 0.0215735 0.965087 600000 0.0200038 0.967603
875000 0.0148509 0.975821 700000 0.0152292 0.975176
1000000 0.00936378 0.984752 800000 0.00896315 0.985384
1125000 0.00490088 0.992007 900000 0.00585213 0.990434
1250000 0.00388501 0.99366 1000000 0.00327599 0.99467
1375000 0.00122887 0.998043 1100000 0.00159845 0.997437
1500000 0.000498063 0.999247 1200000 0.000171108 0.999804
1625000 -0.000468616 1.00085 1300000 -0.000313183 1.00059
1750000 -0.00131731 1.00226 1400000 -0.00148 1.00254
1875000 -0.00247035 1.00419 1500000 -0.00297976 1.00504
2000000 -0.00306799 1.00519 1600000 -0.00287926 1.00487
2125000 -0.00355565 1.006 1700000 -0.00430671 1.00727
2250000 -0.00434289 1.00733 1800000 -0.00453729 1.00765
2375000 -0.00497634 1.00839 1900000 -0.00507356 1.00856
2500000 -0.00576373 1.00972 2000000 -0.00586954 1.0099

53
examples/USER/fep/CH4hyd/fdti01/in.fdti01.lmp
View File

@ -30,51 +30,54 @@ pair_coeff 1 4 lj/cut/tip4p/long/soft 0.0000 1.0000 0.0 # C4H Hw
pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 0.0 # H Ow pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 0.0 # H Ow
pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 0.0 # H Hw pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 0.0 # H Hw
variable nsteps equal 2500000 variable TK equal 300.0
variable nprint equal ${nsteps}/500 variable PBAR equal 1.0
variable ndump equal ${nsteps}/100
variable temp equal 300.0 fix SHAKE all shake 0.0001 20 0 b 2 a 2
variable press equal 1.0
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 a 2
neighbor 2.0 bin neighbor 2.0 bin
timestep 2.0 timestep 1.0
velocity all create ${temp} 12345 velocity all create ${TK} 12345
thermo_style multi thermo_style multi
thermo ${nprint} thermo 5000
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 1000 fix TPSTAT all npt temp ${TK} ${TK} 100 iso ${PBAR} ${PBAR} 1000
run 250000 set type 1*2 charge 0.0
run 100000
reset_timestep 0 reset_timestep 0
variable lambda equal ramp(0.0,1.0) variable lambda equal ramp(0.0,1.0)
variable q1 equal -0.24*v_lambda
variable q2 equal 0.06*v_lambda
fix fADAPT all adapt/fep 125000 pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda after yes fix ADAPT all adapt/fep 100000 &
pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda &
atom charge 1 v_q1 &
atom charge 2 v_q2 &
after yes
fix PRINT all print 100000 "adapt lambda = ${lambda} q1 = ${q1} q2 = ${q2}"
variable dlambda equal 0.002 variable dlambda equal 0.002
variable dq1 equal -0.24*v_dlambda
variable dq2 equal 0.06*v_dlambda
compute cFEP all fep ${temp} pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda compute FEP all fep ${TK} &
pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fdti01.lmp fix FEP all ave/time 20 4000 100000 c_FEP[1] c_FEP[2] file fdti01.lmp
# compute cRDF all rdf 100 3 3 3 4 4 4 dump TRAJ all custom 20000 dump.lammpstrj id mol type element x y z ix iy iz
# fix fRDF all ave/time 2000 500 ${nsteps} c_cRDF file rdf.lammps mode vector dump_modify TRAJ element C H O H
group owater type 3 run 2000000
compute cMSD owater msd
fix fMSD owater ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C H O H
run ${nsteps}
# write_restart restart.*.lmp # write_restart restart.*.lmp
write_data data.*.lmp write_data data.*.lmp

7682
examples/USER/fep/CH4hyd/fdti01/log.lammps
View File

File diff suppressed because it is too large Load Diff

44
examples/USER/fep/CH4hyd/fdti10/fdti10.lmp
View File

@ -1,22 +1,22 @@
# Time-averaged data for fix fFEP # Time-averaged data for fix FEP
# TimeStep c_cFEP[1] c_cFEP[2] # TimeStep c_FEP[1] c_FEP[2]
125000 0.00652155 0.989123 100000 0.00671766 0.988799
250000 0.00572718 0.990445 200000 0.00591179 0.990139
375000 0.00506788 0.991544 300000 0.00495535 0.991733
500000 0.00448303 0.992522 400000 0.00428164 0.992859
625000 0.00370796 0.993818 500000 0.00367253 0.993879
750000 0.00277535 0.995385 600000 0.00334988 0.994424
875000 0.00230833 0.996171 700000 0.00246906 0.995908
1000000 0.000982551 0.998409 800000 0.00139088 0.997725
1125000 0.000314678 0.99954 900000 0.000491084 0.99924
1250000 -0.000839025 1.0015 1000000 0.000615407 0.999044
1375000 -0.00155615 1.00273 1100000 -0.00172714 1.00302
1500000 -0.00320598 1.00553 1200000 -0.00333962 1.00579
1625000 -0.00587133 1.01011 1300000 -0.0040424 1.00699
1750000 -0.00867306 1.01495 1400000 -0.0102385 1.01763
1875000 -0.0134768 1.02327 1500000 -0.0173611 1.03006
2000000 -0.0229314 1.03995 1600000 -0.0243837 1.04234
2125000 -0.0282333 1.04896 1700000 -0.0211518 1.03654
2250000 -0.0202762 1.03477 1800000 -0.0220119 1.03776
2375000 -0.0123279 1.02096 1900000 -0.013304 1.02261
2500000 -0.00546203 1.00923 2000000 -0.00648381 1.01095

52
examples/USER/fep/CH4hyd/fdti10/in.fdti10.lmp
View File

@ -30,51 +30,53 @@ pair_coeff 1 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # C4H Hw
pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 1.0 # H Ow pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 1.0 # H Ow
pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # H Hw pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # H Hw
variable nsteps equal 2500000 variable TK equal 300.0
variable nprint equal ${nsteps}/500 variable PBAR equal 1.0
variable ndump equal ${nsteps}/100
variable temp equal 300.0 fix SHAKE all shake 0.0001 20 0 b 2 a 2
variable press equal 1.0
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 a 2
neighbor 2.0 bin neighbor 2.0 bin
timestep 2.0 timestep 1.0
velocity all create ${temp} 12345 velocity all create ${TK} 12345
thermo_style multi thermo_style multi
thermo ${nprint} thermo 5000
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 1000 fix TPSTAT all npt temp ${TK} ${TK} 100 iso ${PBAR} ${PBAR} 1000
run 250000 run 100000
reset_timestep 0 reset_timestep 0
variable lambda equal ramp(1.0,0.0) variable lambda equal ramp(1.0,0.0)
variable q1 equal -0.24*v_lambda
variable q2 equal 0.06*v_lambda
fix fADAPT all adapt/fep 125000 pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda after yes fix ADAPT all adapt/fep 100000 &
pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda &
atom charge 1 v_q1 &
atom charge 2 v_q2 &
after yes
fix PRINT all print 100000 "adapt lambda = ${lambda} q1 = ${q1} q2 = ${q2}"
variable dlambda equal -0.002 variable dlambda equal -0.002
variable dq1 equal -0.24*v_dlambda
variable dq2 equal 0.06*v_dlambda
compute cFEP all fep ${temp} pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda compute FEP all fep ${TK} &
pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fdti10.lmp fix FEP all ave/time 20 4000 100000 c_FEP[1] c_FEP[2] file fdti10.lmp
# compute cRDF all rdf 100 3 3 3 4 4 4 dump TRAJ all custom 20000 dump.lammpstrj id mol type element x y z ix iy iz
# fix fRDF all ave/time 2000 500 ${nsteps} c_cRDF file rdf.lammps mode vector dump_modify TRAJ element C H O H
group owater type 3 run 2000000
compute cMSD owater msd
fix fMSD owater ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C H O H
run ${nsteps}
# write_restart restart.*.lmp # write_restart restart.*.lmp
write_data data.*.lmp write_data data.*.lmp

7674
examples/USER/fep/CH4hyd/fdti10/log.lammps
View File

File diff suppressed because it is too large Load Diff

44
examples/USER/fep/CH4hyd/fep01/fep01.lmp
View File

@ -1,22 +1,22 @@
# Time-averaged data for fix fFEP # Time-averaged data for fix FEP
# TimeStep c_cFEP[1] c_cFEP[2] # TimeStep c_FEP[1] c_FEP[2]
125000 0.0806308 0.881044 100000 0.0735182 0.889583
250000 0.251798 0.66945 200000 0.241868 0.679931
375000 0.448135 0.502774 300000 0.407677 0.542008
500000 0.670275 0.384034 400000 0.709112 0.360902
625000 0.873864 0.343148 500000 0.718538 0.428553
750000 0.689533 0.484145 600000 0.639674 0.516854
875000 0.474153 0.600305 700000 0.482835 0.586307
1000000 0.303566 0.740249 800000 0.289216 0.746055
1125000 0.165693 0.855143 900000 0.192641 0.823932
1250000 0.130949 0.88524 1000000 0.113029 0.908737
1375000 0.0518948 0.976114 1100000 0.0619301 0.96572
1500000 0.0283335 1.0025 1200000 0.0197356 1.01976
1625000 -0.000790228 1.04259 1300000 0.00310596 1.03223
1750000 -0.0263791 1.0736 1400000 -0.0300484 1.08295
1875000 -0.0587947 1.12568 1500000 -0.0714914 1.14599
2000000 -0.0764999 1.15472 1600000 -0.0712604 1.14573
2125000 -0.091264 1.17869 1700000 -0.109089 1.2123
2250000 -0.11294 1.21778 1800000 -0.117256 1.22671
2375000 -0.130721 1.25162 1900000 -0.132337 1.25534
2500000 -0.151808 1.29392 2000000 -0.153557 1.29745

53
examples/USER/fep/CH4hyd/fep01/in.fep01.lmp
View File

@ -30,51 +30,54 @@ pair_coeff 1 4 lj/cut/tip4p/long/soft 0.0000 1.0000 0.0 # C4H Hw
pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 0.0 # H Ow pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 0.0 # H Ow
pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 0.0 # H Hw pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 0.0 # H Hw
variable nsteps equal 2500000 variable TK equal 300.0
variable nprint equal ${nsteps}/500 variable PBAR equal 1.0
variable ndump equal ${nsteps}/100
variable temp equal 300.0 fix SHAKE all shake 0.0001 20 0 b 2 a 2
variable press equal 1.0
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 a 2
neighbor 2.0 bin neighbor 2.0 bin
timestep 2.0 timestep 1.0
velocity all create ${temp} 12345 velocity all create ${TK} 12345
thermo_style multi thermo_style multi
thermo ${nprint} thermo 5000
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 1000 fix TPSTAT all npt temp ${TK} ${TK} 100 iso ${PBAR} ${PBAR} 1000
run 250000 set type 1*2 charge 0.0
run 100000
reset_timestep 0 reset_timestep 0
variable lambda equal ramp(0.0,1.0) variable lambda equal ramp(0.0,1.0)
variable q1 equal -0.24*v_lambda
variable q2 equal 0.06*v_lambda
fix fADAPT all adapt/fep 125000 pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda after yes fix ADAPT all adapt/fep 100000 &
pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda &
atom charge 1 v_q1 &
atom charge 2 v_q2 &
after yes
fix PRINT all print 100000 "adapt lambda = ${lambda} q1 = ${q1} q2 = ${q2}"
variable dlambda equal 0.05 variable dlambda equal 0.05
variable dq1 equal -0.24*v_dlambda
variable dq2 equal 0.06*v_dlambda
compute cFEP all fep ${temp} pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda compute FEP all fep ${TK} &
pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fep01.lmp fix FEP all ave/time 20 4000 100000 c_FEP[1] c_FEP[2] file fep01.lmp
# compute cRDF all rdf 100 3 3 3 4 4 4 dump TRAJ all custom 20000 dump.lammpstrj id mol type element x y z ix iy iz
# fix fRDF all ave/time 2000 500 ${nsteps} c_cRDF file rdf.lammps mode vector dump_modify TRAJ element C H O H
group owater type 3 run 2000000
compute cMSD owater msd
fix fMSD owater ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C H O H
run ${nsteps}
# write_restart restart.*.lmp # write_restart restart.*.lmp
write_data data.*.lmp write_data data.*.lmp

7682
examples/USER/fep/CH4hyd/fep01/log.lammps
View File

File diff suppressed because it is too large Load Diff

44
examples/USER/fep/CH4hyd/fep10/fep10.lmp
View File

@ -1,22 +1,22 @@
# Time-averaged data for fix fFEP # Time-averaged data for fix FEP
# TimeStep c_cFEP[1] c_cFEP[2] # TimeStep c_FEP[1] c_FEP[2]
125000 0.156203 0.771285 100000 0.160455 0.766272
250000 0.136448 0.799233 200000 0.141912 0.791902
375000 0.12198 0.820426 300000 0.119585 0.824744
500000 0.109453 0.840255 400000 0.104694 0.847809
625000 0.0921141 0.868102 500000 0.0917124 0.869318
750000 0.0713015 0.905735 600000 0.0859541 0.881974
875000 0.0622431 0.922352 700000 0.0666319 0.920268
1000000 0.0330546 0.981942 800000 0.0428654 0.965893
1125000 0.0197866 1.01134 900000 0.0240106 0.999889
1250000 -0.00412826 1.06765 1000000 0.0294147 0.997671
1375000 -0.0174302 1.11343 1100000 -0.0214947 1.11583
1500000 -0.0511646 1.18453 1200000 -0.0543642 1.23078
1625000 -0.105676 1.36325 1300000 -0.0665772 1.25986
1750000 -0.161616 1.56478 1400000 -0.196675 1.66521
1875000 -0.258263 1.93825 1500000 -0.341037 2.39926
2000000 -0.447487 3.0595 1600000 -0.480217 2.97287
2125000 -0.549377 3.01163 1700000 -0.405599 2.36311
2250000 -0.379009 2.01933 1800000 -0.415165 2.12414
2375000 -0.20934 1.45678 1900000 -0.229669 1.49644
2500000 -0.0644056 1.1252 2000000 -0.0838847 1.15743

52
examples/USER/fep/CH4hyd/fep10/in.fep10.lmp
View File

@ -30,51 +30,53 @@ pair_coeff 1 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # C4H Hw
pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 1.0 # H Ow pair_coeff 2 3 lj/cut/tip4p/long/soft 0.0699 2.8126 1.0 # H Ow
pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # H Hw pair_coeff 2 4 lj/cut/tip4p/long/soft 0.0000 1.0000 1.0 # H Hw
variable nsteps equal 2500000 variable TK equal 300.0
variable nprint equal ${nsteps}/500 variable PBAR equal 1.0
variable ndump equal ${nsteps}/100
variable temp equal 300.0 fix SHAKE all shake 0.0001 20 0 b 2 a 2
variable press equal 1.0
fix fSHAKE all shake 0.0001 20 ${nprint} b 2 a 2
neighbor 2.0 bin neighbor 2.0 bin
timestep 2.0 timestep 1.0
velocity all create ${temp} 12345 velocity all create ${TK} 12345
thermo_style multi thermo_style multi
thermo ${nprint} thermo 5000
fix fNPT all npt temp ${temp} ${temp} 100 iso ${press} ${press} 1000 fix TPSTAT all npt temp ${TK} ${TK} 100 iso ${PBAR} ${PBAR} 1000
run 250000 run 100000
reset_timestep 0 reset_timestep 0
variable lambda equal ramp(1.0,0.0) variable lambda equal ramp(1.0,0.0)
variable q1 equal -0.24*v_lambda
variable q2 equal 0.06*v_lambda
fix fADAPT all adapt/fep 125000 pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda after yes fix ADAPT all adapt/fep 100000 &
pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_lambda &
atom charge 1 v_q1 &
atom charge 2 v_q2 &
after yes
fix PRINT all print 100000 "adapt lambda = ${lambda} q1 = ${q1} q2 = ${q2}"
variable dlambda equal -0.05 variable dlambda equal -0.05
variable dq1 equal -0.24*v_dlambda
variable dq2 equal 0.06*v_dlambda
compute cFEP all fep ${temp} pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda compute FEP all fep ${TK} &
pair lj/cut/tip4p/long/soft lambda 1*2 3*4 v_dlambda &
atom charge 1 v_dq1 &
atom charge 2 v_dq2
fix fFEP all ave/time 25 4000 125000 c_cFEP[1] c_cFEP[2] file fep10.lmp fix FEP all ave/time 20 4000 100000 c_FEP[1] c_FEP[2] file fep10.lmp
# compute cRDF all rdf 100 3 3 3 4 4 4 dump TRAJ all custom 20000 dump.lammpstrj id mol type element x y z ix iy iz
# fix fRDF all ave/time 2000 500 ${nsteps} c_cRDF file rdf.lammps mode vector dump_modify TRAJ element C H O H
group owater type 3 run 2000000
compute cMSD owater msd
fix fMSD owater ave/time 1 1 ${ndump} c_cMSD[1] c_cMSD[2] c_cMSD[3] c_cMSD[4] file msd.lammps
dump dCONF all custom ${ndump} dump.lammpstrj id mol type element x y z ix iy iz
dump_modify dCONF element C H O H
run ${nsteps}
# write_restart restart.*.lmp # write_restart restart.*.lmp
write_data data.*.lmp write_data data.*.lmp

7674
examples/USER/fep/CH4hyd/fep10/log.lammps
View File

File diff suppressed because it is too large Load Diff

24
lib/gpu/lal_coul_dsf.cu
View File

@ -81,7 +81,7 @@ __kernel void k_coul_dsf(const __global numtyp4 *restrict x_,
int j=dev_packed[nbor]; int j=dev_packed[nbor];
numtyp factor_coul, r, prefactor, erfcc; numtyp factor_coul, r, prefactor, erfcc;
factor_coul = sp_lj[sbmask(j)]; factor_coul = (numtyp)1.0-sp_lj[sbmask(j)];
j &= NEIGHMASK; j &= NEIGHMASK;
numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j]; numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
@ -98,12 +98,12 @@ __kernel void k_coul_dsf(const __global numtyp4 *restrict x_,
r = ucl_sqrt(rsq); r = ucl_sqrt(rsq);
fetch(prefactor,j,q_tex); fetch(prefactor,j,q_tex);
prefactor *= factor_coul * qqrd2e*qtmp/r; prefactor *= qqrd2e*qtmp/r;
numtyp erfcd = ucl_exp(-alpha*alpha*rsq); numtyp erfcd = ucl_exp(-alpha*alpha*rsq);
numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*alpha*r); numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc + (numtyp)2.0*alpha/MY_PIS*r*erfcd + forcecoul = prefactor * (erfcc + (numtyp)2.0*alpha/MY_PIS*r*erfcd +
rsq*f_shift); rsq*f_shift-factor_coul);
force = forcecoul * r2inv; force = forcecoul * r2inv;
@ -112,10 +112,8 @@ __kernel void k_coul_dsf(const __global numtyp4 *restrict x_,
f.z+=delz*force; f.z+=delz*force;
if (eflag>0) { if (eflag>0) {
if (rsq < cut_coulsq) { numtyp e=prefactor*(erfcc-r*e_shift-rsq*f_shift-factor_coul);
numtyp e=prefactor*(erfcc-r*e_shift-rsq*f_shift); e_coul += e;
e_coul += e;
}
} }
if (vflag>0) { if (vflag>0) {
virial[0] += delx*delx*force; virial[0] += delx*delx*force;
@ -182,7 +180,7 @@ __kernel void k_coul_dsf_fast(const __global numtyp4 *restrict x_,
int j=dev_packed[nbor]; int j=dev_packed[nbor];
numtyp factor_coul, r, prefactor, erfcc; numtyp factor_coul, r, prefactor, erfcc;
factor_coul = sp_lj[sbmask(j)]; factor_coul = (numtyp)1.0-sp_lj[sbmask(j)];
j &= NEIGHMASK; j &= NEIGHMASK;
numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j]; numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
@ -199,12 +197,12 @@ __kernel void k_coul_dsf_fast(const __global numtyp4 *restrict x_,
r = ucl_sqrt(rsq); r = ucl_sqrt(rsq);
fetch(prefactor,j,q_tex); fetch(prefactor,j,q_tex);
prefactor *= factor_coul * qqrd2e*qtmp/r; prefactor *= qqrd2e*qtmp/r;
numtyp erfcd = ucl_exp(-alpha*alpha*rsq); numtyp erfcd = ucl_exp(-alpha*alpha*rsq);
numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*alpha*r); numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc + (numtyp)2.0*alpha/MY_PIS*r*erfcd + forcecoul = prefactor * (erfcc + (numtyp)2.0*alpha/MY_PIS*r*erfcd +
rsq*f_shift); rsq*f_shift-factor_coul);
force = forcecoul * r2inv; force = forcecoul * r2inv;
@ -213,10 +211,8 @@ __kernel void k_coul_dsf_fast(const __global numtyp4 *restrict x_,
f.z+=delz*force; f.z+=delz*force;
if (eflag>0) { if (eflag>0) {
if (rsq < cut_coulsq) { numtyp e=prefactor*(erfcc-r*e_shift-rsq*f_shift-factor_coul);
numtyp e=prefactor*(erfcc-r*e_shift-rsq*f_shift); e_coul += e;
e_coul += e;
}
} }
if (vflag>0) { if (vflag>0) {
virial[0] += delx*delx*force; virial[0] += delx*delx*force;

16
lib/gpu/lal_lj_dsf.cu
View File

@ -89,7 +89,7 @@ __kernel void k_lj_dsf(const __global numtyp4 *restrict x_,
numtyp factor_lj, factor_coul, r, prefactor, erfcc; numtyp factor_lj, factor_coul, r, prefactor, erfcc;
factor_lj = sp_lj[sbmask(j)]; factor_lj = sp_lj[sbmask(j)];
factor_coul = sp_lj[sbmask(j)+4]; factor_coul = (numtyp)1.0-sp_lj[sbmask(j)+4];
j &= NEIGHMASK; j &= NEIGHMASK;
numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j]; numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
@ -115,12 +115,12 @@ __kernel void k_lj_dsf(const __global numtyp4 *restrict x_,
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
r = ucl_sqrt(rsq); r = ucl_sqrt(rsq);
fetch(prefactor,j,q_tex); fetch(prefactor,j,q_tex);
prefactor *= factor_coul * qqrd2e*qtmp/r; prefactor *= qqrd2e*qtmp/r;
numtyp erfcd = ucl_exp(-alpha*alpha*rsq); numtyp erfcd = ucl_exp(-alpha*alpha*rsq);
numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*alpha*r); numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc + (numtyp)2.0*alpha/MY_PIS*r*erfcd + forcecoul = prefactor * (erfcc + (numtyp)2.0*alpha/MY_PIS*r*erfcd +
rsq*f_shift); rsq*f_shift-factor_coul);
} else } else
forcecoul = (numtyp)0.0; forcecoul = (numtyp)0.0;
@ -132,7 +132,7 @@ __kernel void k_lj_dsf(const __global numtyp4 *restrict x_,
if (eflag>0) { if (eflag>0) {
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
numtyp e=prefactor*(erfcc-r*e_shift-rsq*f_shift); numtyp e=prefactor*(erfcc-r*e_shift-rsq*f_shift-factor_coul);
e_coul += e; e_coul += e;
} }
if (rsq < lj1[mtype].z) { if (rsq < lj1[mtype].z) {
@ -217,7 +217,7 @@ __kernel void k_lj_dsf_fast(const __global numtyp4 *restrict x_,
numtyp factor_lj, factor_coul, r, prefactor, erfcc; numtyp factor_lj, factor_coul, r, prefactor, erfcc;
factor_lj = sp_lj[sbmask(j)]; factor_lj = sp_lj[sbmask(j)];
factor_coul = sp_lj[sbmask(j)+4]; factor_coul = (numtyp)1.0-sp_lj[sbmask(j)+4];
j &= NEIGHMASK; j &= NEIGHMASK;
numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j]; numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
@ -242,12 +242,12 @@ __kernel void k_lj_dsf_fast(const __global numtyp4 *restrict x_,
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
r = ucl_sqrt(rsq); r = ucl_sqrt(rsq);
fetch(prefactor,j,q_tex); fetch(prefactor,j,q_tex);
prefactor *= factor_coul * qqrd2e*qtmp/r; prefactor *= qqrd2e*qtmp/r;
numtyp erfcd = ucl_exp(-alpha*alpha*rsq); numtyp erfcd = ucl_exp(-alpha*alpha*rsq);
numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*alpha*r); numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc + (numtyp)2.0*alpha/MY_PIS*r*erfcd + forcecoul = prefactor * (erfcc + (numtyp)2.0*alpha/MY_PIS*r*erfcd +
rsq*f_shift); rsq*f_shift-factor_coul);
} else } else
forcecoul = (numtyp)0.0; forcecoul = (numtyp)0.0;
@ -259,7 +259,7 @@ __kernel void k_lj_dsf_fast(const __global numtyp4 *restrict x_,
if (eflag>0) { if (eflag>0) {
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
numtyp e=prefactor*(erfcc-r*e_shift-rsq*f_shift); numtyp e=prefactor*(erfcc-r*e_shift-rsq*f_shift-factor_coul);
e_coul += e; e_coul += e;
} }
if (rsq < lj1[mtype].z) { if (rsq < lj1[mtype].z) {

10
src/GPU/pair_coul_dsf_gpu.cpp
View File

@ -203,7 +203,7 @@ void PairCoulDSFGPU::cpu_compute(int start, int inum, int eflag, int vflag,
int i,j,ii,jj,jnum; int i,j,ii,jj,jnum;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
double r,rsq,r2inv,forcecoul,factor_coul; double r,rsq,r2inv,forcecoul,factor_coul;
double prefactor,erfcc,erfcd,e_self,t; double prefactor,erfcc,erfcd,t;
int *jlist; int *jlist;
ecoul = 0.0; ecoul = 0.0;
@ -227,7 +227,7 @@ void PairCoulDSFGPU::cpu_compute(int start, int inum, int eflag, int vflag,
jnum = numneigh[i]; jnum = numneigh[i];
if (evflag) { if (evflag) {
e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e; double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e;
ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0); ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0);
} }
@ -244,13 +244,14 @@ void PairCoulDSFGPU::cpu_compute(int start, int inum, int eflag, int vflag,
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
r2inv = 1.0/rsq; r2inv = 1.0/rsq;
r = sqrt(rsq); r = sqrt(rsq);
prefactor = factor_coul * qqrd2e*qtmp*q[j]/r; prefactor = qqrd2e*qtmp*q[j]/r;
erfcd = exp(-alpha*alpha*r*r); erfcd = exp(-alpha*alpha*r*r);
t = 1.0 / (1.0 + EWALD_P*alpha*r); t = 1.0 / (1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd + forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
r*f_shift) * r; r*f_shift) * r;
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
fpair = forcecoul * r2inv; fpair = forcecoul * r2inv;
f[i][0] += delx*fpair; f[i][0] += delx*fpair;
f[i][1] += dely*fpair; f[i][1] += dely*fpair;
@ -259,6 +260,7 @@ void PairCoulDSFGPU::cpu_compute(int start, int inum, int eflag, int vflag,
if (eflag) { if (eflag) {
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0; } else ecoul = 0.0;
} }

8
src/GPU/pair_lj_cut_coul_dsf_gpu.cpp
View File

@ -205,7 +205,7 @@ void PairLJCutCoulDSFGPU::cpu_compute(int start, int inum, int eflag, int vflag,
int i,j,ii,jj,jnum,itype,jtype; int i,j,ii,jj,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
double prefactor,erfcc,erfcd,e_self,t; double prefactor,erfcc,erfcd,t;
int *jlist; int *jlist;
evdwl = ecoul = 0.0; evdwl = ecoul = 0.0;
@ -232,7 +232,7 @@ void PairLJCutCoulDSFGPU::cpu_compute(int start, int inum, int eflag, int vflag,
jnum = numneigh[i]; jnum = numneigh[i];
if (evflag) { if (evflag) {
e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e; double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e;
ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0); ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0);
} }
@ -258,12 +258,13 @@ void PairLJCutCoulDSFGPU::cpu_compute(int start, int inum, int eflag, int vflag,
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
r = sqrt(rsq); r = sqrt(rsq);
prefactor = factor_coul * qqrd2e*qtmp*q[j]/r; prefactor = qqrd2e*qtmp*q[j]/r;
erfcd = exp(-alpha*alpha*r*r); erfcd = exp(-alpha*alpha*r*r);
t = 1.0 / (1.0 + EWALD_P*alpha*r); t = 1.0 / (1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd + forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
r*f_shift) * r; r*f_shift) * r;
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} }
fpair = (forcecoul + factor_lj*forcelj) * r2inv; fpair = (forcecoul + factor_lj*forcelj) * r2inv;
@ -280,6 +281,7 @@ void PairLJCutCoulDSFGPU::cpu_compute(int start, int inum, int eflag, int vflag,
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0; } else ecoul = 0.0;
} }

2
src/USER-OMP/pair_born_coul_wolf_omp.cpp
View File

@ -125,7 +125,7 @@ void PairBornCoulWolfOMP::eval(int iifrom, int iito, ThrData * const thr)
qisq = qtmp*qtmp; qisq = qtmp*qtmp;
e_self = -(e_shift/2.0 + alf/MY_PIS) * qisq*qqrd2e; e_self = -(e_shift/2.0 + alf/MY_PIS) * qisq*qqrd2e;
if (EVFLAG) ev_tally_thr(this,i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0,thr); if (EFLAG) ev_tally_thr(this,i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0,thr);
for (jj = 0; jj < jnum; jj++) { for (jj = 0; jj < jnum; jj++) {
j = jlist[jj]; j = jlist[jj];

9
src/USER-OMP/pair_coul_dsf_omp.cpp
View File

@ -138,14 +138,14 @@ void PairCoulDSFOMP::eval(int iifrom, int iito, ThrData * const thr)
r2inv = 1.0/rsq; r2inv = 1.0/rsq;
r = sqrt(rsq); r = sqrt(rsq);
prefactor = factor_coul * qqrd2e*qtmp*q[j]/r; prefactor = qqrd2e*qtmp*q[j]/r;
erfcd = exp(-alpha*alpha*rsq); erfcd = exp(-alpha*alpha*rsq);
t = 1.0 / (1.0 + EWALD_P*alpha*r); t = 1.0 / (1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd + forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
r*f_shift) * r; r*f_shift) * r;
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
fpair = forcecoul * r2inv; fpair = forcecoul * r2inv;
if (EFLAG) ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
fxtmp += delx*fpair; fxtmp += delx*fpair;
fytmp += dely*fpair; fytmp += dely*fpair;
@ -156,6 +156,11 @@ void PairCoulDSFOMP::eval(int iifrom, int iito, ThrData * const thr)
f[j].z -= delz*fpair; f[j].z -= delz*fpair;
} }
if (EFLAG) {
ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR, if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
0.0,ecoul,fpair,delx,dely,delz,thr); 0.0,ecoul,fpair,delx,dely,delz,thr);
} }

21
src/USER-OMP/pair_lj_cut_coul_dsf_omp.cpp
View File

@ -122,7 +122,7 @@ void PairLJCutCoulDSFOMP::eval(int iifrom, int iito, ThrData * const thr)
jnum = numneigh[i]; jnum = numneigh[i];
fxtmp=fytmp=fztmp=0.0; fxtmp=fytmp=fztmp=0.0;
if (EVFLAG) { if (EFLAG) {
double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e; double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e;
ev_tally_thr(this,i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0,thr); ev_tally_thr(this,i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0,thr);
} }
@ -145,19 +145,21 @@ void PairLJCutCoulDSFOMP::eval(int iifrom, int iito, ThrData * const thr)
if (rsq < cut_ljsq[itype][jtype]) { if (rsq < cut_ljsq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv; r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
forcelj *= factor_lj;
} else forcelj = 0.0; } else forcelj = 0.0;
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
r = sqrt(rsq); r = sqrt(rsq);
prefactor = factor_coul * qqrd2e*qtmp*q[j]/r; prefactor = qqrd2e*qtmp*q[j]/r;
erfcd = exp(-alpha*alpha*r*r); erfcd = exp(-alpha*alpha*r*r);
t = 1.0 / (1.0 + EWALD_P*alpha*r); t = 1.0 / (1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd + forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
r*f_shift) * r; r*f_shift) * r;
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
fpair = forcecoul * r2inv;
} else forcecoul = 0.0; } else forcecoul = 0.0;
fpair = (forcecoul + forcelj) * r2inv;
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
fxtmp += delx*fpair; fxtmp += delx*fpair;
fytmp += dely*fpair; fytmp += dely*fpair;
@ -171,12 +173,13 @@ void PairLJCutCoulDSFOMP::eval(int iifrom, int iito, ThrData * const thr)
if (EFLAG) { if (EFLAG) {
if (rsq < cut_ljsq[itype][jtype]) { if (rsq < cut_ljsq[itype][jtype]) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype]; offset[itype][jtype];
evdwl *= factor_lj; evdwl *= factor_lj;
} else evdwl = 0.0; } else evdwl = 0.0;
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0; } else ecoul = 0.0;
} }

2
src/pair_born_coul_wolf.cpp
View File

@ -115,7 +115,7 @@ void PairBornCoulWolf::compute(int eflag, int vflag)
qisq = qtmp*qtmp; qisq = qtmp*qtmp;
e_self = -(e_shift/2.0 + alf/MY_PIS) * qisq*qqrd2e; e_self = -(e_shift/2.0 + alf/MY_PIS) * qisq*qqrd2e;
if (evflag) ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0); if (eflag) ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0);
for (jj = 0; jj < jnum; jj++) { for (jj = 0; jj < jnum; jj++) {
j = jlist[jj]; j = jlist[jj];

10
src/pair_coul_dsf.cpp
View File

@ -67,7 +67,7 @@ void PairCoulDSF::compute(int eflag, int vflag)
int i,j,ii,jj,inum,jnum; int i,j,ii,jj,inum,jnum;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
double r,rsq,r2inv,forcecoul,factor_coul; double r,rsq,r2inv,forcecoul,factor_coul;
double prefactor,erfcc,erfcd,e_self,t; double prefactor,erfcc,erfcd,t;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
ecoul = 0.0; ecoul = 0.0;
@ -99,7 +99,7 @@ void PairCoulDSF::compute(int eflag, int vflag)
jnum = numneigh[i]; jnum = numneigh[i];
if (eflag) { if (eflag) {
e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e; double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e;
ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0); ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0);
} }
@ -117,14 +117,15 @@ void PairCoulDSF::compute(int eflag, int vflag)
r2inv = 1.0/rsq; r2inv = 1.0/rsq;
r = sqrt(rsq); r = sqrt(rsq);
prefactor = factor_coul * qqrd2e*qtmp*q[j]/r; prefactor = qqrd2e*qtmp*q[j]/r;
erfcd = exp(-alpha*alpha*rsq); erfcd = exp(-alpha*alpha*rsq);
t = 1.0 / (1.0 + EWALD_P*alpha*r); t = 1.0 / (1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd + forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
r*f_shift) * r; r*f_shift) * r;
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
fpair = forcecoul * r2inv; fpair = forcecoul * r2inv;
f[i][0] += delx*fpair; f[i][0] += delx*fpair;
f[i][1] += dely*fpair; f[i][1] += dely*fpair;
f[i][2] += delz*fpair; f[i][2] += delz*fpair;
@ -136,6 +137,7 @@ void PairCoulDSF::compute(int eflag, int vflag)
if (eflag) { if (eflag) {
ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0; } else ecoul = 0.0;
if (evflag) ev_tally(i,j,nlocal,newton_pair, if (evflag) ev_tally(i,j,nlocal,newton_pair,

8
src/pair_coul_wolf.cpp
View File

@ -133,11 +133,9 @@ void PairCoulWolf::compute(int eflag, int vflag)
} }
if (eflag) { if (eflag) {
if (rsq < cut_coulsq) { ecoul = v_sh;
ecoul = v_sh; if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; } else ecoul = 0.0;
} else ecoul = 0.0;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair, if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,ecoul,fpair,delx,dely,delz); 0.0,ecoul,fpair,delx,dely,delz);

9
src/pair_lj_cut_coul_dsf.cpp
View File

@ -138,12 +138,14 @@ void PairLJCutCoulDSF::compute(int eflag, int vflag)
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
r = sqrt(rsq); r = sqrt(rsq);
prefactor = factor_coul * qqrd2e*qtmp*q[j]/r; prefactor = qqrd2e*qtmp*q[j]/r;
erfcd = exp(-alpha*alpha*r*r); erfcd = exp(-alpha*alpha*r*r);
t = 1.0 / (1.0 + EWALD_P*alpha*r); t = 1.0 / (1.0 + EWALD_P*alpha*r);
erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd + forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
r*f_shift) * r; r*f_shift) * r;
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
fpair = forcecoul * r2inv;
} else forcecoul = 0.0; } else forcecoul = 0.0;
fpair = (forcecoul + factor_lj*forcelj) * r2inv; fpair = (forcecoul + factor_lj*forcelj) * r2inv;
@ -162,9 +164,10 @@ void PairLJCutCoulDSF::compute(int eflag, int vflag)
offset[itype][jtype]; offset[itype][jtype];
evdwl *= factor_lj; evdwl *= factor_lj;
} else evdwl = 0.0; } else evdwl = 0.0;
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0; } else ecoul = 0.0;
} }