diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index a88fa6f6bf..bfa07041b3 100644 --- a/doc/src/bond_bpm_rotational.rst +++ b/doc/src/bond_bpm_rotational.rst @@ -68,16 +68,14 @@ Bonds will break under sufficient stress. A breaking criteria is calculated .. math:: - B = \alpha(r, r_0) \frac{f_r}{f_{r,c}} + \frac{|f_s|}{f_{s,c}} + - \frac{|\tau_b|}{\tau_{b,c}} + \frac{|\tau_t|}{\tau_{t,c}} + B = \mathrm{max}\{0, \frac{f_r}{f_{r,c}} + \frac{|f_s|}{f_{s,c}} + + \frac{|\tau_b|}{\tau_{b,c}} + \frac{|\tau_t|}{\tau_{t,c}} \} where :math:`|f_s|` is the magnitude of the shear force and :math:`|\tau_b|` and :math:`|\tau_t|` are the magnitudes of the bending and twisting forces, respectively. The corresponding variables :math:`f_{r,c}` :math:`f_{s,c}`, :math:`\tau_{b,c}`, and :math:`\tau_{t,c}` are critical -limits to each force or torque. The term :math:`\alpha` is simply one in -extension and zero in compression such that the normal force component -does not contribute to the breaking criteria in compression. +limits to each force or torque. If :math:`B` is ever equal to or exceeds one, the bond will break. This is done by setting by setting its type to 0 such that forces and torques are no longer computed. diff --git a/doc/src/read_data.rst b/doc/src/read_data.rst index 5b5c951688..e4cf2eff0f 100644 --- a/doc/src/read_data.rst +++ b/doc/src/read_data.rst @@ -648,6 +648,8 @@ of analysis. - atom-ID atom-type rho esph cv x y z * - sphere - atom-ID atom-type diameter density x y z + * - sphere/bpm + - atom-ID molecule-ID atom-type diameter density x y z * - spin - atom-ID atom-type x y z spx spy spz sp * - tdpd diff --git a/examples/bpm/impact/impact_bpm.lmp b/examples/bpm/impact/impact_bpm.lmp index 9eb04829bb..de975fb956 100644 --- a/examples/bpm/impact/impact_bpm.lmp +++ b/examples/bpm/impact/impact_bpm.lmp @@ -23,8 +23,8 @@ neighbor 1.0 bin pair_style gran/hooke/history 1.0 NULL 0.5 NULL 0.1 1 bond_style bpm/rotational store/local 2 time id1 id2 pair_coeff 1 1 -bond_coeff 1 1.0 0.2 0.02 0.02 0.05 0.01 0.005 0.005 0.1 0.02 0.002 0.002 -bond_coeff 2 1.0 0.2 0.02 0.02 0.20 0.04 0.020 0.020 0.1 0.02 0.002 0.002 +bond_coeff 1 1.0 0.2 0.02 0.02 0.05 0.01 0.01 0.01 0.1 0.02 0.002 0.002 +bond_coeff 2 1.0 0.2 0.02 0.02 0.20 0.04 0.04 0.04 0.1 0.02 0.002 0.002 fix 1 all nve/sphere/bpm fix 2 all store/local 100 3 @@ -49,4 +49,4 @@ thermo 100 dump 2 all local 100 brokenDump f_2[1] f_2[2] f_2[3] dump_modify 2 header no -run 6000 +run 10000 diff --git a/examples/bpm/impact/log.30Jul2021.pour.g++1.4 b/examples/bpm/impact/log.30Jul2021.pour.g++1.4 index 50d90788c6..300ac271eb 100644 --- a/examples/bpm/impact/log.30Jul2021.pour.g++1.4 +++ b/examples/bpm/impact/log.30Jul2021.pour.g++1.4 @@ -36,8 +36,8 @@ neighbor 1.0 bin pair_style gran/hooke/history 1.0 NULL 0.5 NULL 0.1 1 bond_style bpm/rotational store/local 2 time id1 id2 pair_coeff 1 1 -bond_coeff 1 1.0 0.2 0.02 0.02 0.05 0.01 0.005 0.005 0.1 0.02 0.002 0.002 -bond_coeff 2 1.0 0.2 0.02 0.02 0.20 0.04 0.020 0.020 0.1 0.02 0.002 0.002 +bond_coeff 1 1.0 0.2 0.02 0.02 0.05 0.01 0.01 0.01 0.1 0.02 0.002 0.002 +bond_coeff 2 1.0 0.2 0.02 0.02 0.20 0.04 0.04 0.04 0.1 0.02 0.002 0.002 fix 1 all nve/sphere/bpm fix 2 all store/local 100 3 @@ -67,7 +67,7 @@ Finding 1-2 1-3 1-4 neighbors ... special bond factors coul: 0 1 1 15 = max # of 1-2 neighbors 101 = max # of special neighbors - special bonds CPU = 0.001 seconds + special bonds CPU = 0.002 seconds create_bonds many projectile projectile 2 0.0 1.5 Added 21869 bonds, new total = 60420 Finding 1-2 1-3 1-4 neighbors ... @@ -94,7 +94,7 @@ dump 1 all custom 100 atomDump id radius x y z c_nbond dump 2 all local 100 brokenDump f_2[1] f_2[2] f_2[3] dump_modify 2 header no -run 6000 +run 10000 Neighbor list info ... update every 1 steps, delay 10 steps, check yes max neighbors/atom: 2000, page size: 100000 @@ -114,89 +114,128 @@ Step KinEng PotEng Pxx Pyy Pzz c_tbond 200 0.00053248439 -5.4354084 0.00013279102 1.6787387e-12 5.4953548e-13 10.870817 300 0.00053248439 -5.4354084 0.00013459514 3.0577302e-12 -1.2088316e-12 10.870817 400 0.00053248439 -5.4354084 0.00013739515 9.7326931e-13 1.1562543e-12 10.870817 - 500 0.0005311195 -5.4353511 0.00013532545 -1.8468866e-06 -1.6585907e-06 10.870817 - 600 0.00051833989 -5.4288362 7.1725873e-05 -1.6444182e-05 -1.9715516e-05 10.860022 - 700 0.00050552886 -5.3985398 -2.3212715e-05 -1.3937832e-05 -1.9091149e-05 10.805955 - 800 0.00050518224 -5.3616112 8.2515921e-05 2.6973287e-05 2.0513255e-05 10.735516 - 900 0.00050396335 -5.3458995 0.00019263348 -3.7564641e-05 -4.5902247e-05 10.699172 - 1000 0.00049399705 -5.3356803 8.3780211e-05 -3.900811e-05 -3.7974569e-05 10.677852 - 1100 0.00047949984 -5.3177445 -3.3039707e-05 7.3302691e-05 9.1425255e-05 10.645196 - 1200 0.0004778277 -5.2918722 2.5525662e-05 -7.3957522e-06 -1.0428816e-06 10.596168 - 1300 0.00047175855 -5.2622229 0.00016622071 -7.5895427e-05 -7.3668463e-05 10.536884 - 1400 0.00047346655 -5.2357068 0.00013088251 3.8629925e-05 4.7370005e-05 10.484257 - 1500 0.00047030222 -5.2110028 2.0138342e-05 0.00014344934 0.00014572537 10.43289 - 1600 0.00047120889 -5.1944302 8.0308725e-05 9.7652665e-05 9.2786939e-05 10.397265 - 1700 0.00046169345 -5.1800484 0.00015698427 -1.2415417e-05 -2.2822524e-05 10.370457 - 1800 0.00045838981 -5.1581008 6.4666951e-05 -2.044593e-05 -3.1699628e-05 10.327726 - 1900 0.00045723383 -5.1435986 7.4830972e-06 2.8528283e-05 9.5685497e-06 10.29615 - 2000 0.00045822947 -5.1339093 8.8355461e-05 2.4391936e-05 2.6856807e-06 10.276268 - 2100 0.00045635704 -5.123239 0.00014447351 -1.885824e-05 -2.5038419e-05 10.252789 - 2200 0.0004531016 -5.1193556 4.9062089e-05 -4.1290134e-05 -2.8125502e-05 10.245682 - 2300 0.00044750986 -5.1059583 -1.3106261e-05 -6.3125644e-05 -3.3287217e-05 10.221303 - 2400 0.00044828631 -5.08466 9.1580692e-05 -3.8583903e-05 -5.1373805e-05 10.177942 - 2500 0.00044674619 -5.0736124 0.00011856396 4.6266072e-06 -1.5907979e-05 10.153203 - 2600 0.00044727275 -5.0614096 5.295308e-05 6.0578952e-05 5.7167864e-05 10.128553 - 2700 0.00044110051 -5.0507956 1.8440804e-05 5.4883922e-05 6.3982178e-05 10.106153 - 2800 0.00043633042 -5.0426937 8.1816485e-05 -1.5431256e-06 -1.855006e-05 10.08996 - 2900 0.00043176823 -5.0356307 0.00012046457 4.4651701e-05 2.7620283e-06 10.075387 - 3000 0.00043218951 -5.027846 3.2574549e-05 7.4433525e-05 4.6834269e-05 10.059554 - 3100 0.00043256828 -5.0161104 3.3101063e-05 4.6551359e-05 5.6531167e-05 10.036164 - 3200 0.0004332929 -5.0068366 8.9512897e-05 -5.7060048e-06 -2.5891142e-06 10.018172 - 3300 0.00043157298 -5.0008275 8.4139723e-05 -3.6940688e-05 -5.1910847e-05 10.005218 - 3400 0.000431398 -4.9965161 2.5625084e-05 1.1433716e-05 6.4926012e-06 9.9965815 - 3500 0.00043011473 -4.9928073 2.7095736e-05 -7.7339585e-06 7.844574e-06 9.9891148 - 3600 0.00042888594 -4.9847529 7.7453653e-05 -4.6984476e-05 -3.5038838e-05 9.9747211 - 3700 0.00042753057 -4.9765311 6.8775078e-05 -3.1077337e-05 -3.3804378e-05 9.9585282 - 3800 0.00042500059 -4.9705379 1.4686991e-05 -2.1799517e-05 -1.2214392e-05 9.9471033 - 3900 0.00042443518 -4.9635965 2.9590319e-05 -7.5174838e-06 9.4138716e-06 9.9336092 - 4000 0.00042354431 -4.9575067 7.497612e-05 -2.241558e-05 -1.7795282e-05 9.9216445 - 4100 0.00042603308 -4.937097 4.8886365e-05 -1.065802e-05 -2.4723829e-05 9.8824217 - 4200 0.00042476916 -4.926145 1.9200413e-05 1.0927038e-05 3.1855248e-06 9.8600216 - 4300 0.00042485728 -4.918146 3.0086129e-05 -8.6292641e-06 2.1921266e-06 9.8440986 - 4400 0.00042494366 -4.9099091 6.2777705e-05 -5.4921818e-06 -7.6169646e-06 9.8263764 - 4500 0.0004275409 -4.8995178 4.7617152e-05 2.1335231e-06 1.1579137e-06 9.8055056 - 4600 0.00042662895 -4.8820208 1.084188e-05 -1.5696264e-06 1.1310734e-05 9.7706909 - 4700 0.00042734689 -4.8725198 3.7425941e-05 -4.439827e-06 1.13299e-05 9.7511695 - 4800 0.00042900821 -4.8540514 5.4751797e-05 -1.1417464e-05 -8.8084486e-06 9.7147355 - 4900 0.00043023872 -4.8363433 3.3053643e-05 3.2468406e-06 -6.4727375e-06 9.6776718 - 5000 0.00043199875 -4.8198092 1.8345516e-05 8.4946117e-06 5.6309681e-06 9.6456459 - 5100 0.00044223594 -4.7778052 3.0243074e-05 -5.0316681e-06 -3.9518237e-06 9.5634221 - 5200 0.00044479909 -4.7640293 4.8598154e-05 -6.7947105e-06 -1.5950295e-05 9.5343649 - 5300 0.0004454718 -4.755521 2.7021591e-05 4.8540854e-06 -5.1116404e-06 9.5168226 - 5400 0.00044509902 -4.7442744 1.5098441e-05 9.1872595e-06 1.0269456e-05 9.4952321 - 5500 0.00044706883 -4.7308865 3.2631779e-05 1.097946e-06 5.8901266e-06 9.4690536 - 5600 0.00045084112 -4.7115078 3.7164777e-05 -2.2594948e-06 -5.1676349e-06 9.4311803 - 5700 0.00045327828 -4.6989042 2.2566508e-05 2.2741586e-06 8.6663706e-07 9.404552 - 5800 0.00046103079 -4.6780119 1.4582664e-05 1.2169813e-06 8.2506998e-06 9.3626304 - 5900 0.00046315911 -4.6539134 2.7960455e-05 -5.7750919e-06 1.1483488e-08 9.313692 - 6000 0.00046381574 -4.6376345 3.2506838e-05 -1.1407228e-07 -5.8332064e-06 9.2804066 -Loop time of 26.8573 on 4 procs for 6000 steps with 11116 atoms + 500 0.00053111833 -5.4353812 0.00013532154 -1.8602155e-06 -1.6702711e-06 10.870817 + 600 0.00051558884 -5.4335908 5.8899364e-05 -2.4651157e-05 -2.6584555e-05 10.868838 + 700 0.0005013055 -5.4216397 -7.255462e-05 -2.8970044e-05 -4.1925574e-05 10.850846 + 800 0.00049878634 -5.3947446 5.6196335e-05 1.8649508e-05 8.9976433e-06 10.803706 + 900 0.00049589843 -5.3619596 0.00020821095 -7.7726675e-05 -8.6749853e-05 10.737855 + 1000 0.00049487709 -5.3376233 0.00017638288 -2.8676423e-05 -4.5487748e-05 10.683879 + 1100 0.00048720924 -5.3283688 4.0171402e-05 0.0001406456 0.00014483 10.662648 + 1200 0.00048274004 -5.3195134 7.4824656e-06 7.4411513e-05 0.00010365246 10.643397 + 1300 0.00047990668 -5.3127593 0.00018114159 -2.9205108e-05 1.1366149e-05 10.629723 + 1400 0.00048034109 -5.3084824 0.000147425 -6.2432251e-05 -3.8277687e-05 10.620187 + 1500 0.00047225373 -5.3024141 8.6252599e-07 2.4082822e-05 3.049023e-05 10.608492 + 1600 0.00045698513 -5.2972292 1.0763242e-05 4.1148987e-05 2.83019e-05 10.599676 + 1700 0.00044854655 -5.2896066 0.00010217014 -6.9430229e-05 -9.6661607e-05 10.587262 + 1800 0.00044929529 -5.2775304 0.00010981512 -8.0553726e-05 -0.00011931185 10.562972 + 1900 0.00044441992 -5.2650251 1.8882842e-05 1.0491309e-05 -3.879472e-05 10.540842 + 2000 0.00043947002 -5.2484561 4.4869915e-05 6.8824059e-05 3.5808833e-05 10.511155 + 2100 0.00043372382 -5.2265524 0.00013970367 1.6484426e-05 2.2785645e-05 10.469054 + 2200 0.00043174768 -5.2069377 9.4662371e-05 2.2278647e-06 3.5962708e-05 10.427852 + 2300 0.0004270123 -5.1924656 2.3188212e-05 6.1522399e-06 7.4262622e-05 10.396006 + 2400 0.00042569191 -5.1848751 5.7811979e-05 2.7124197e-05 0.0001072424 10.378733 + 2500 0.00042259717 -5.1758217 0.00013069307 3.3656662e-05 0.00010197635 10.359482 + 2600 0.00042271375 -5.1700793 9.5183077e-05 5.3232578e-05 8.9912385e-05 10.347787 + 2700 0.00042252395 -5.1628759 2.6305437e-05 6.7644983e-05 7.7095778e-05 10.332674 + 2800 0.00042123755 -5.1549973 6.5875753e-05 2.8392569e-05 2.8079356e-05 10.315941 + 2900 0.00042132346 -5.1508594 0.00011264272 6.3410829e-06 -2.1983564e-05 10.306405 + 3000 0.00042343054 -5.1495182 7.3503612e-05 3.8773748e-05 -9.519275e-06 10.301547 + 3100 0.00042220412 -5.1483831 2.5233575e-05 5.0076476e-05 3.938746e-06 10.299568 + 3200 0.00042303815 -5.1466902 6.1939664e-05 -9.1346169e-06 -2.5920911e-05 10.29651 + 3300 0.00042000178 -5.144782 9.8814555e-05 -4.04749e-05 -4.2876825e-05 10.292371 + 3400 0.00041874209 -5.1416065 5.1246647e-05 -2.7877246e-05 -3.225052e-05 10.286254 + 3500 0.00041582277 -5.1397016 2.0053694e-05 -3.5797833e-05 -1.536015e-05 10.282296 + 3600 0.00041607097 -5.139236 6.0675623e-05 -5.7232123e-05 -3.1162791e-05 10.281036 + 3700 0.00041445536 -5.1373913 8.7909083e-05 -4.1136114e-05 -5.2627526e-05 10.277978 + 3800 0.0004147342 -5.1323122 4.6324048e-05 7.0253754e-06 -3.3511914e-05 10.268442 + 3900 0.00041446917 -5.1294358 2.8646507e-05 1.5201733e-05 -1.13759e-05 10.262504 + 4000 0.00041346205 -5.1250314 6.540586e-05 -1.17595e-05 -2.8050171e-05 10.254948 + 4100 0.00041230785 -5.1219436 7.9364924e-05 -4.1504333e-06 -2.2530525e-05 10.248111 + 4200 0.00041198555 -5.1177883 4.3515184e-05 1.5227343e-05 -6.3707934e-06 10.240014 + 4300 0.0004111489 -5.1134893 2.8350236e-05 2.0718016e-06 1.2010375e-05 10.231558 + 4400 0.00041090623 -5.1104369 6.2460747e-05 -2.5959985e-05 7.8242641e-07 10.224901 + 4500 0.00040944466 -5.1085221 6.7135567e-05 -1.7699087e-05 -4.7022089e-06 10.220943 + 4600 0.00040810594 -5.1065034 3.5212952e-05 1.3568365e-05 1.3875486e-05 10.217704 + 4700 0.00040810646 -5.1039941 3.3409499e-05 1.7215022e-05 2.8204859e-05 10.212487 + 4800 0.0004074837 -5.1015741 5.5792503e-05 1.9563116e-06 1.823506e-06 10.207089 + 4900 0.00040677077 -5.0987121 5.6695901e-05 8.2729584e-06 -1.2713008e-05 10.201871 + 5000 0.00040636045 -5.0961728 3.0704198e-05 2.2141861e-05 8.2099332e-06 10.196474 + 5100 0.00040606831 -5.0947673 3.1281394e-05 7.0864149e-06 2.0262936e-05 10.193055 + 5200 0.00040652265 -5.0940213 5.2610835e-05 -1.2888854e-05 3.0894446e-06 10.191076 + 5300 0.00040642029 -5.0931407 4.6148958e-05 -9.5544284e-06 -6.047443e-06 10.189457 + 5400 0.00040642806 -5.0915733 2.59528e-05 -3.3035524e-06 1.026995e-05 10.186038 + 5500 0.00040686546 -5.0908 2.9026708e-05 -9.0382082e-06 1.4643294e-05 10.184059 + 5600 0.0004064361 -5.0908057 4.6731327e-05 -1.2664731e-05 -2.6232887e-06 10.183339 + 5700 0.00040629203 -5.0903044 4.0856223e-05 -1.2201759e-06 -1.3169401e-05 10.18262 + 5800 0.00040637688 -5.0890571 2.2625249e-05 8.7645385e-06 -6.2486963e-06 10.180281 + 5900 0.00040613721 -5.0874767 2.8382883e-05 2.4072343e-06 2.0419388e-07 10.176862 + 6000 0.00040668084 -5.0865465 4.3602089e-05 -5.4962058e-06 -4.5087421e-06 10.174523 + 6100 0.00040707325 -5.0865389 3.4958823e-05 -2.0750124e-06 -1.6708517e-06 10.174343 + 6200 0.00040691768 -5.0863974 2.1602821e-05 3.1566836e-06 1.0526645e-05 10.174343 + 6300 0.00040705673 -5.0862887 2.9412395e-05 9.2283412e-07 1.4273225e-05 10.173983 + 6400 0.00040648035 -5.0860509 3.9684464e-05 -1.663237e-06 3.9771927e-06 10.173624 + 6500 0.00040710623 -5.0861041 3.1078617e-05 1.732822e-06 3.4003409e-06 10.173624 + 6600 0.00040665879 -5.0857907 2.0693771e-05 3.3053846e-06 1.2181329e-05 10.173264 + 6700 0.00040650151 -5.0854203 2.8479998e-05 2.7244033e-06 9.7566436e-06 10.172364 + 6800 0.00040635649 -5.0851816 3.5150661e-05 1.6906684e-06 -6.155957e-06 10.171644 + 6900 0.00040670804 -5.0848974 2.4632227e-05 2.9367001e-06 -1.0691056e-05 10.171105 + 7000 0.00040693354 -5.0843039 1.8740741e-05 6.024808e-06 4.2065619e-07 10.169665 + 7100 0.00040728228 -5.0843393 2.7137965e-05 5.5506702e-06 5.5908974e-06 10.169665 + 7200 0.0004074084 -5.0842875 3.0307934e-05 2.0170793e-06 -2.4296651e-06 10.169485 + 7300 0.00040723509 -5.0843468 2.1465618e-05 1.9642493e-06 -3.6022271e-06 10.169485 + 7400 0.00040756027 -5.0843623 1.6801323e-05 -1.9748948e-06 1.4205345e-06 10.169306 + 7500 0.00040829979 -5.0843202 2.4772881e-05 -6.1251363e-06 2.0247483e-06 10.169126 + 7600 0.00040822994 -5.0843182 2.7272667e-05 -3.0357928e-06 4.894101e-07 10.169126 + 7700 0.00040831723 -5.0843052 1.9410405e-05 3.6094291e-06 1.5451381e-06 10.169126 + 7800 0.00040868149 -5.0843706 1.6484224e-05 3.3901782e-06 3.9911363e-06 10.169126 + 7900 0.00040872521 -5.0843735 2.2844838e-05 2.8813595e-06 1.4869802e-06 10.169126 + 8000 0.00040853749 -5.0843239 2.3537039e-05 5.1951501e-06 -1.2448734e-06 10.169126 + 8100 0.00040812899 -5.0842554 1.6947117e-05 7.5128919e-06 -1.0877933e-06 10.168946 + 8200 0.00040812313 -5.0842813 1.5639254e-05 3.6719094e-06 -9.3497783e-07 10.168946 + 8300 0.00040817027 -5.0842752 2.0634335e-05 2.5358492e-07 -3.2726598e-06 10.168946 + 8400 0.00040774138 -5.084215 2.0224447e-05 1.3696075e-06 -3.3568279e-06 10.168766 + 8500 0.00040760502 -5.0842 1.4541525e-05 9.3556598e-07 1.1823477e-06 10.168766 + 8600 0.00040756971 -5.0841463 1.4460781e-05 -2.7822738e-06 4.3070092e-06 10.168766 + 8700 0.00040706312 -5.0840255 1.8278276e-05 -5.20189e-06 1.0784628e-06 10.168766 + 8800 0.00040670111 -5.0839094 1.7116511e-05 -9.4769204e-07 -3.2089738e-06 10.168586 + 8900 0.00040617439 -5.0838164 1.3315166e-05 3.2313582e-06 -2.3144297e-06 10.168586 + 9000 0.00040576758 -5.0837468 1.3898828e-05 1.5947021e-06 -2.0719014e-06 10.168586 + 9100 0.00040577217 -5.0837244 1.6547097e-05 1.1667189e-06 -3.2056138e-06 10.168406 + 9200 0.00040615545 -5.0837984 1.4946269e-05 4.3601683e-06 -2.1585248e-06 10.168406 + 9300 0.00040638526 -5.083836 1.1737091e-05 5.1607613e-06 7.2161152e-07 10.168406 + 9400 0.0004062125 -5.0838558 1.2486756e-05 1.9996225e-06 1.6192477e-06 10.168406 + 9500 0.00040627984 -5.0839239 1.441806e-05 -6.6274154e-07 -2.9396969e-07 10.168406 + 9600 0.0004065461 -5.0839109 1.3189089e-05 -5.1486848e-07 4.6653236e-07 10.168406 + 9700 0.00040642188 -5.0838722 1.0626956e-05 -1.7580535e-06 2.8200944e-06 10.168226 + 9800 0.0004061705 -5.0838326 1.1280802e-05 -3.4868266e-06 2.7287279e-06 10.168226 + 9900 0.00040666798 -5.0835647 1.2432396e-05 -2.8727154e-06 1.4556152e-07 10.167686 + 10000 0.00040675506 -5.0831833 1.0832242e-05 4.3061564e-07 -4.1422229e-07 10.166967 +Loop time of 48.2068 on 4 procs for 10000 steps with 11116 atoms -Performance: 965099.761 tau/day, 223.403 timesteps/s -98.4% CPU use with 4 MPI tasks x no OpenMP threads +Performance: 896139.501 tau/day, 207.440 timesteps/s +97.7% CPU use with 4 MPI tasks x no OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 0.24862 | 0.26655 | 0.28715 | 3.1 | 0.99 -Bond | 18.756 | 20.055 | 21.432 | 25.3 | 74.67 -Neigh | 0.68191 | 0.68601 | 0.69029 | 0.4 | 2.55 -Comm | 1.6547 | 3.1706 | 4.5875 | 69.7 | 11.81 -Output | 0.3381 | 0.33879 | 0.34058 | 0.2 | 1.26 -Modify | 2.0658 | 2.1974 | 2.3512 | 7.6 | 8.18 -Other | | 0.1428 | | | 0.53 +Pair | 0.39096 | 0.42232 | 0.45519 | 4.6 | 0.88 +Bond | 30.522 | 33.954 | 38.166 | 51.7 | 70.43 +Neigh | 1.1822 | 1.1872 | 1.1915 | 0.3 | 2.46 +Comm | 3.4201 | 7.9545 | 11.664 | 115.6 | 16.50 +Output | 1.0078 | 1.0085 | 1.0099 | 0.1 | 2.09 +Modify | 3.0972 | 3.3841 | 3.7307 | 14.3 | 7.02 +Other | | 0.2958 | | | 0.61 -Nlocal: 2779.00 ave 3444 max 2189 min -Histogram: 1 0 1 0 0 0 1 0 0 1 -Nghost: 1160.25 ave 1335 max 985 min -Histogram: 1 0 1 0 0 0 0 1 0 1 -Neighs: 10709.5 ave 13715 max 8243 min +Nlocal: 2779.00 ave 4159 max 1635 min Histogram: 2 0 0 0 0 0 0 1 0 1 +Nghost: 1002.00 ave 1229 max 800 min +Histogram: 1 1 0 0 0 0 0 1 0 1 +Neighs: 11431.8 ave 18381 max 6059 min +Histogram: 2 0 0 0 0 0 1 0 0 1 -Total # of neighbors = 42838 -Ave neighs/atom = 3.8537244 -Ave special neighs/atom = 9.2786974 -Neighbor list builds = 402 +Total # of neighbors = 45727 +Ave neighs/atom = 4.1136200 +Ave special neighs/atom = 10.166967 +Neighbor list builds = 637 Dangerous builds = 0 - -Total wall time: 0:00:26 +Total wall time: 0:00:48 diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index 7a0d4fedd6..a6e5d4583c 100644 --- a/src/BPM/bond_bpm.cpp +++ b/src/BPM/bond_bpm.cpp @@ -1,349 +1,346 @@ -// clang-format off -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "bond_bpm.h" - -#include "atom.h" -#include "domain.h" -#include "error.h" -#include "fix_store_local.h" -#include "fix_update_special_bonds.h" -#include "force.h" -#include "memory.h" -#include "modify.h" -#include "update.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp) -{ - id_fix_store_local = nullptr; - id_fix_prop_atom = nullptr; - fix_store_local = nullptr; - fix_update_special_bonds = nullptr; - - prop_atom_flag = 0; - nvalues = 0; - output_data = nullptr; - pack_choice = nullptr; - - r0_max_estimate = 0.0; - max_stretch = 1.0; -} - -/* ---------------------------------------------------------------------- */ - -BondBPM::~BondBPM() -{ - delete [] pack_choice; - delete [] id_fix_store_local; - delete [] id_fix_prop_atom; - memory->destroy(output_data); -} - -/* ---------------------------------------------------------------------- */ - -void BondBPM::init_style() -{ - int ifix; - if (id_fix_store_local) { - ifix = modify->find_fix(id_fix_store_local); - if (ifix < 0) error->all(FLERR, "Cannot find fix store/local"); - if (strcmp(modify->fix[ifix]->style, "store/local") != 0) - error->all(FLERR, "Incorrect fix style matched, not store/local"); - fix_store_local = (FixStoreLocal *) modify->fix[ifix]; - fix_store_local->nvalues = nvalues; - } - - ifix = modify->find_fix_by_style("update/special/bonds"); - if (ifix >= 0) - fix_update_special_bonds = (FixUpdateSpecialBonds *) modify->fix[ifix]; - else - fix_update_special_bonds = nullptr; - - - if (force->angle || force->dihedral || force->improper) - error->all(FLERR, - "Bond style bpm cannot be used with 3,4-body interactions"); - if (atom->molecular == 2) - error->all(FLERR, - "Bond style bpm cannot be used with atom style template"); - - // special 1-3 and 1-4 weights must be 1 to prevent building 1-3 and 1-4 special bond lists - if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || - force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) - error->all(FLERR,"Bond style bpm requires 1-3 and 1-4 special weights of 1.0"); -} - -/* ---------------------------------------------------------------------- - global settings - All args before store/local command are saved for potential args - for specific bond BPM substyles - All args after optional store/local command are variables stored - in the compute store/local -------------------------------------------------------------------------- */ - -void BondBPM::settings(int narg, char **arg) -{ - int iarg = 0; - while (iarg < narg) { - if (id_fix_store_local) { - if (strcmp(arg[iarg], "id1") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_id1; - } else if (strcmp(arg[iarg], "id2") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_id2; - } else if (strcmp(arg[iarg], "time") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_time; - } else if (strcmp(arg[iarg], "x") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_x; - } else if (strcmp(arg[iarg], "y") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_y; - } else if (strcmp(arg[iarg], "z") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_z; - } else if (strcmp(arg[iarg], "x/ref") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_x_ref; - prop_atom_flag = 1; - } else if (strcmp(arg[iarg], "y/ref") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_y_ref; - prop_atom_flag = 1; - } else if (strcmp(arg[iarg], "z/ref") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_z_ref; - prop_atom_flag = 1; - } else { - error->all(FLERR, "Illegal bond_style command"); - } - } else if (strcmp(arg[iarg], "store/local") == 0) { - id_fix_store_local = utils::strdup(arg[iarg+1]); - iarg ++; - nvalues = 0; - pack_choice = new FnPtrPack[narg - iarg - 1]; - } - iarg ++; - } - - if (id_fix_store_local) { - if (nvalues == 0) error->all(FLERR, - "Bond style bpm/rotational must include at least one value to output"); - memory->create(output_data, nvalues, "bond/bpm:output_data"); - - // Use store property to save reference positions as it can transfer to ghost atoms - if (prop_atom_flag == 1) { - - id_fix_prop_atom = utils::strdup("BPM_PROPERTY_ATOM"); - int ifix = modify->find_fix(id_fix_prop_atom); - if (ifix < 0) { - modify->add_fix(fmt::format("{} all property/atom " - "d_BPM_X_REF d_BPM_Y_REF d_BPM_Z_REF ghost yes", id_fix_prop_atom)); - ifix = modify->find_fix(id_fix_prop_atom); - } - - int type_flag; - int col_flag; - index_x_ref = atom->find_custom("BPM_X_REF", type_flag, col_flag); - index_y_ref = atom->find_custom("BPM_Y_REF", type_flag, col_flag); - index_z_ref = atom->find_custom("BPM_Z_REF", type_flag, col_flag); - - if (modify->fix[ifix]->restart_reset) { - modify->fix[ifix]->restart_reset = 0; - } else { - double *x_ref = atom->dvector[index_x_ref]; - double *y_ref = atom->dvector[index_y_ref]; - double *z_ref = atom->dvector[index_z_ref]; - - double **x = atom->x; - for (int i = 0; i < atom->nlocal; i++) { - x_ref[i] = x[i][0]; - y_ref[i] = x[i][1]; - z_ref[i] = x[i][2]; - } - } - } - } -} - -/* ---------------------------------------------------------------------- - used to check bond communiction cutoff - not perfect, estimates based on local-local only -------------------------------------------------------------------------- */ - -double BondBPM::equilibrium_distance(int i) -{ - // Ghost atoms may not yet be communicated, this may only be an estimate - if (r0_max_estimate == 0) { - int type, j; - double delx, dely, delz, r; - double **x = atom->x; - for (int i = 0; i < atom->nlocal; i ++) { - for (int m = 0; m < atom->num_bond[i]; m ++) { - type = atom->bond_type[i][m]; - if (type == 0) continue; - - j = atom->map(atom->bond_atom[i][m]); - if(j == -1) continue; - - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - domain->minimum_image(delx, dely, delz); - - r = sqrt(delx*delx + dely*dely + delz*delz); - if(r > r0_max_estimate) r0_max_estimate = r; - } - } - - double temp; - MPI_Allreduce(&r0_max_estimate,&temp,1,MPI_DOUBLE,MPI_MAX,world); - r0_max_estimate = temp; - - //if (comm->me == 0) - // utils::logmesg(lmp,fmt::format("Estimating longest bond = {}\n",r0_max_estimate)); - } - - // Divide out heuristic prefactor added in comm class - return max_stretch*r0_max_estimate/1.5; -} - -/* ---------------------------------------------------------------------- */ - -void BondBPM::process_broken(int i, int j) -{ - if (fix_store_local) { - for (int n = 0; n < nvalues; n++) - (this->*pack_choice[n])(n, i, j); - - fix_store_local->add_data(output_data, i, j); - } - - if (fix_update_special_bonds) - fix_update_special_bonds->add_broken_bond(i, j); - - // Manually search and remove from atom arrays - // need to remove in case special bonds arrays rebuilt - int m, n; - int nlocal = atom->nlocal; - - tagint *tag = atom->tag; - tagint **bond_atom = atom->bond_atom; - int **bond_type = atom->bond_type; - int *num_bond = atom->num_bond; - - if (i < nlocal) { - for (m = 0; m < num_bond[i]; m++) { - // if(tag[i] == 25 or tag[j] == 25) printf("\t 1 ATOM LOOP %d-%d, %d/%d \n", tag[i], bond_atom[i][m], m, num_bond[i]); - if (bond_atom[i][m] == tag[j]) { - // if(tag[i] == 25 or tag[j] == 25) printf("\t DELETED\n"); - bond_type[i][m] = 0; - n = num_bond[i]; - bond_type[i][m] = bond_type[i][n-1]; - bond_atom[i][m] = bond_atom[i][n-1]; - num_bond[i]--; - break; - } - } - } - - if (j < nlocal) { - for (m = 0; m < num_bond[j]; m++) { - // if(tag[i] == 25 or tag[j] == 25) printf("\t 2 ATOM LOOP %d-%d, %d/%d \n", tag[j], bond_atom[j][m], m, num_bond[j]); - if (bond_atom[j][m] == tag[i]) { - // if(tag[j] == 25 or tag[j] == 25) printf("\t DELETED\n"); - bond_type[j][m] = 0; - n = num_bond[j]; - bond_type[j][m] = bond_type[j][n-1]; - bond_atom[j][m] = bond_atom[j][n-1]; - num_bond[j]--; - break; - } - } - } - - // if((tag[i] == 10 and tag[j] == 23)or (tag[i] == 23 and tag[j] == 10)) printf("Broken bond %d-%d (%d %d/%d)\n", tag[i], tag[j], i, j, nlocal); - -} - -/* ---------------------------------------------------------------------- - one method for every keyword bond bpm can output - the atom property is packed into array or vector -------------------------------------------------------------------------- */ - -void BondBPM::pack_id1(int n, int i, int j) -{ - tagint *tag = atom->tag; - output_data[n] = tag[i]; -} - -/* ---------------------------------------------------------------------- */ - -void BondBPM::pack_id2(int n, int i, int j) -{ - tagint *tag = atom->tag; - output_data[n] = tag[j]; -} - -/* ---------------------------------------------------------------------- */ - -void BondBPM::pack_time(int n, int i, int j) -{ - bigint time = update->ntimestep; - output_data[n] = time; -} - -/* ---------------------------------------------------------------------- */ - -void BondBPM::pack_x(int n, int i, int j) -{ - double **x = atom->x; - output_data[n] = (x[i][0] + x[j][0])*0.5; -} - -/* ---------------------------------------------------------------------- */ - -void BondBPM::pack_y(int n, int i, int j) -{ - double **x = atom->x; - output_data[n] = (x[i][1] + x[j][1])*0.5; -} - -/* ---------------------------------------------------------------------- */ - -void BondBPM::pack_z(int n, int i, int j) -{ - double **x = atom->x; - output_data[n] = (x[i][2] + x[j][2])*0.5; -} - -/* ---------------------------------------------------------------------- */ - -void BondBPM::pack_x_ref(int n, int i, int j) -{ - double *x = atom->dvector[index_x_ref]; - output_data[n] = (x[i] + x[j])*0.5; -} - -/* ---------------------------------------------------------------------- */ - -void BondBPM::pack_y_ref(int n, int i, int j) -{ - double *y = atom->dvector[index_y_ref]; - output_data[n] = (y[i] + y[j])*0.5; -} - -/* ---------------------------------------------------------------------- */ - -void BondBPM::pack_z_ref(int n, int i, int j) -{ - double *z = atom->dvector[index_z_ref]; - output_data[n] = (z[i] + z[j])*0.5; -} +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "bond_bpm.h" + +#include "atom.h" +#include "domain.h" +#include "error.h" +#include "fix_bond_history.h" +#include "fix_store_local.h" +#include "fix_update_special_bonds.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp) +{ + id_fix_store_local = nullptr; + id_fix_prop_atom = nullptr; + fix_store_local = nullptr; + fix_update_special_bonds = nullptr; + fix_bond_history = nullptr; + + prop_atom_flag = 0; + nvalues = 0; + output_data = nullptr; + pack_choice = nullptr; + + r0_max_estimate = 0.0; + max_stretch = 1.0; +} + +/* ---------------------------------------------------------------------- */ + +BondBPM::~BondBPM() +{ + delete [] pack_choice; + delete [] id_fix_store_local; + delete [] id_fix_prop_atom; + memory->destroy(output_data); +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::init_style() +{ + int ifix; + if (id_fix_store_local) { + ifix = modify->find_fix(id_fix_store_local); + if (ifix < 0) error->all(FLERR, "Cannot find fix store/local"); + if (strcmp(modify->fix[ifix]->style, "store/local") != 0) + error->all(FLERR, "Incorrect fix style matched, not store/local"); + fix_store_local = (FixStoreLocal *) modify->fix[ifix]; + fix_store_local->nvalues = nvalues; + } + + ifix = modify->find_fix_by_style("update/special/bonds"); + if (ifix >= 0) + fix_update_special_bonds = (FixUpdateSpecialBonds *) modify->fix[ifix]; + else + fix_update_special_bonds = nullptr; + + + if (force->angle || force->dihedral || force->improper) + error->all(FLERR, + "Bond style bpm cannot be used with 3,4-body interactions"); + if (atom->molecular == 2) + error->all(FLERR, + "Bond style bpm cannot be used with atom style template"); + + // special 1-3 and 1-4 weights must be 1 to prevent building 1-3 and 1-4 special bond lists + if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || + force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) + error->all(FLERR,"Bond style bpm requires 1-3 and 1-4 special weights of 1.0"); +} + +/* ---------------------------------------------------------------------- + global settings + All args before store/local command are saved for potential args + for specific bond BPM substyles + All args after optional store/local command are variables stored + in the compute store/local +------------------------------------------------------------------------- */ + +void BondBPM::settings(int narg, char **arg) +{ + int iarg = 0; + while (iarg < narg) { + if (id_fix_store_local) { + if (strcmp(arg[iarg], "id1") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_id1; + } else if (strcmp(arg[iarg], "id2") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_id2; + } else if (strcmp(arg[iarg], "time") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_time; + } else if (strcmp(arg[iarg], "x") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_x; + } else if (strcmp(arg[iarg], "y") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_y; + } else if (strcmp(arg[iarg], "z") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_z; + } else if (strcmp(arg[iarg], "x/ref") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_x_ref; + prop_atom_flag = 1; + } else if (strcmp(arg[iarg], "y/ref") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_y_ref; + prop_atom_flag = 1; + } else if (strcmp(arg[iarg], "z/ref") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_z_ref; + prop_atom_flag = 1; + } else { + error->all(FLERR, "Illegal bond_style command"); + } + } else if (strcmp(arg[iarg], "store/local") == 0) { + id_fix_store_local = utils::strdup(arg[iarg+1]); + iarg ++; + nvalues = 0; + pack_choice = new FnPtrPack[narg - iarg - 1]; + } + iarg ++; + } + + if (id_fix_store_local) { + if (nvalues == 0) error->all(FLERR, + "Bond style bpm/rotational must include at least one value to output"); + memory->create(output_data, nvalues, "bond/bpm:output_data"); + + // Use store property to save reference positions as it can transfer to ghost atoms + if (prop_atom_flag == 1) { + + id_fix_prop_atom = utils::strdup("BPM_PROPERTY_ATOM"); + int ifix = modify->find_fix(id_fix_prop_atom); + if (ifix < 0) { + modify->add_fix(fmt::format("{} all property/atom " + "d_BPM_X_REF d_BPM_Y_REF d_BPM_Z_REF ghost yes", id_fix_prop_atom)); + ifix = modify->find_fix(id_fix_prop_atom); + } + + int type_flag; + int col_flag; + index_x_ref = atom->find_custom("BPM_X_REF", type_flag, col_flag); + index_y_ref = atom->find_custom("BPM_Y_REF", type_flag, col_flag); + index_z_ref = atom->find_custom("BPM_Z_REF", type_flag, col_flag); + + if (modify->fix[ifix]->restart_reset) { + modify->fix[ifix]->restart_reset = 0; + } else { + double *x_ref = atom->dvector[index_x_ref]; + double *y_ref = atom->dvector[index_y_ref]; + double *z_ref = atom->dvector[index_z_ref]; + + double **x = atom->x; + for (int i = 0; i < atom->nlocal; i++) { + x_ref[i] = x[i][0]; + y_ref[i] = x[i][1]; + z_ref[i] = x[i][2]; + } + } + } + } +} + +/* ---------------------------------------------------------------------- + used to check bond communiction cutoff - not perfect, estimates based on local-local only +------------------------------------------------------------------------- */ + +double BondBPM::equilibrium_distance(int i) +{ + // Ghost atoms may not yet be communicated, this may only be an estimate + if (r0_max_estimate == 0) { + int type, j; + double delx, dely, delz, r; + double **x = atom->x; + for (int i = 0; i < atom->nlocal; i ++) { + for (int m = 0; m < atom->num_bond[i]; m ++) { + type = atom->bond_type[i][m]; + if (type == 0) continue; + + j = atom->map(atom->bond_atom[i][m]); + if(j == -1) continue; + + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + domain->minimum_image(delx, dely, delz); + + r = sqrt(delx*delx + dely*dely + delz*delz); + if(r > r0_max_estimate) r0_max_estimate = r; + } + } + + double temp; + MPI_Allreduce(&r0_max_estimate,&temp,1,MPI_DOUBLE,MPI_MAX,world); + r0_max_estimate = temp; + + //if (comm->me == 0) + // utils::logmesg(lmp,fmt::format("Estimating longest bond = {}\n",r0_max_estimate)); + } + + // Divide out heuristic prefactor added in comm class + return max_stretch*r0_max_estimate/1.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::process_broken(int i, int j) +{ + if (fix_store_local) { + for (int n = 0; n < nvalues; n++) + (this->*pack_choice[n])(n, i, j); + + fix_store_local->add_data(output_data, i, j); + } + + if (fix_update_special_bonds) + fix_update_special_bonds->add_broken_bond(i, j); + + // Manually search and remove from atom arrays + // need to remove in case special bonds arrays rebuilt + int m, n; + int nlocal = atom->nlocal; + + tagint *tag = atom->tag; + tagint **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + int *num_bond = atom->num_bond; + + if (i < nlocal) { + for (m = 0; m < num_bond[i]; m++) { + if (bond_atom[i][m] == tag[j]) { + bond_type[i][m] = 0; + n = num_bond[i]; + bond_type[i][m] = bond_type[i][n-1]; + bond_atom[i][m] = bond_atom[i][n-1]; + fix_bond_history->delete_bond(i, m); + num_bond[i]--; + break; + } + } + } + + if (j < nlocal) { + for (m = 0; m < num_bond[j]; m++) { + if (bond_atom[j][m] == tag[i]) { + bond_type[j][m] = 0; + n = num_bond[j]; + bond_type[j][m] = bond_type[j][n-1]; + bond_atom[j][m] = bond_atom[j][n-1]; + fix_bond_history->delete_bond(j, m); + num_bond[j]--; + break; + } + } + } +} + +/* ---------------------------------------------------------------------- + one method for every keyword bond bpm can output + the atom property is packed into array or vector +------------------------------------------------------------------------- */ + +void BondBPM::pack_id1(int n, int i, int j) +{ + tagint *tag = atom->tag; + output_data[n] = tag[i]; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_id2(int n, int i, int j) +{ + tagint *tag = atom->tag; + output_data[n] = tag[j]; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_time(int n, int i, int j) +{ + bigint time = update->ntimestep; + output_data[n] = time; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_x(int n, int i, int j) +{ + double **x = atom->x; + output_data[n] = (x[i][0] + x[j][0])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_y(int n, int i, int j) +{ + double **x = atom->x; + output_data[n] = (x[i][1] + x[j][1])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_z(int n, int i, int j) +{ + double **x = atom->x; + output_data[n] = (x[i][2] + x[j][2])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_x_ref(int n, int i, int j) +{ + double *x = atom->dvector[index_x_ref]; + output_data[n] = (x[i] + x[j])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_y_ref(int n, int i, int j) +{ + double *y = atom->dvector[index_y_ref]; + output_data[n] = (y[i] + y[j])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_z_ref(int n, int i, int j) +{ + double *z = atom->dvector[index_z_ref]; + output_data[n] = (z[i] + z[j])*0.5; +} diff --git a/src/BPM/bond_bpm.h b/src/BPM/bond_bpm.h index d77e44e1d6..53a202802d 100644 --- a/src/BPM/bond_bpm.h +++ b/src/BPM/bond_bpm.h @@ -1,90 +1,91 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef LMP_BOND_BPM_H -#define LMP_BOND_BPM_H - -#include "bond.h" - -namespace LAMMPS_NS { - -class BondBPM : public Bond { - public: - BondBPM(class LAMMPS *); - virtual ~BondBPM(); - virtual void compute(int, int) = 0; - virtual void coeff(int, char **) = 0; - virtual void init_style(); - void settings(int, char **); - double equilibrium_distance(int); - void write_restart(FILE *){}; - void read_restart(FILE *){}; - void write_data(FILE *) {}; - double single(int, double, int, int, double &) = 0; - - protected: - double r0_max_estimate; - double max_stretch; - - char *id_fix_store_local, *id_fix_prop_atom; - class FixStoreLocal *fix_store_local; - class FixUpdateSpecialBonds *fix_update_special_bonds; - - void process_broken(int, int); - typedef void (BondBPM::*FnPtrPack)(int,int,int); - FnPtrPack *pack_choice; // ptrs to pack functions - double *output_data; - - int prop_atom_flag, nvalues; - int index_x_ref, index_y_ref, index_z_ref; - - void pack_id1(int,int,int); - void pack_id2(int,int,int); - void pack_time(int,int,int); - void pack_x(int,int,int); - void pack_y(int,int,int); - void pack_z(int,int,int); - void pack_x_ref(int,int,int); - void pack_y_ref(int,int,int); - void pack_z_ref(int,int,int); -}; - -} // namespace LAMMPS_NS - -#endif - -/* ERROR/WARNING messages: - -E: Cannot find fix store/local - -Fix id cannot be found. - -E: Illegal bond_style command - -Self-explanatory. - -E: Bond style bpm/rotational must include at least one value to output - -Must include at least one bond property to store in fix store/local - -E: Bond style bpm/rotational cannot be used with 3,4-body interactions - -No angle, dihedral, or improper styles can be defined when using -bond style bpm/rotational. - -E: Bond style bpm/rotational cannot be used with atom style template - -This bond style can change the bond topology which is not -allowed with this atom style. - -*/ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_BOND_BPM_H +#define LMP_BOND_BPM_H + +#include "bond.h" + +namespace LAMMPS_NS { + +class BondBPM : public Bond { + public: + BondBPM(class LAMMPS *); + virtual ~BondBPM(); + virtual void compute(int, int) = 0; + virtual void coeff(int, char **) = 0; + virtual void init_style(); + void settings(int, char **); + double equilibrium_distance(int); + void write_restart(FILE *){}; + void read_restart(FILE *){}; + void write_data(FILE *) {}; + double single(int, double, int, int, double &) = 0; + + protected: + double r0_max_estimate; + double max_stretch; + + char *id_fix_store_local, *id_fix_prop_atom; + class FixStoreLocal *fix_store_local; + class FixUpdateSpecialBonds *fix_update_special_bonds; + class FixBondHistory *fix_bond_history; + + void process_broken(int, int); + typedef void (BondBPM::*FnPtrPack)(int,int,int); + FnPtrPack *pack_choice; // ptrs to pack functions + double *output_data; + + int prop_atom_flag, nvalues; + int index_x_ref, index_y_ref, index_z_ref; + + void pack_id1(int,int,int); + void pack_id2(int,int,int); + void pack_time(int,int,int); + void pack_x(int,int,int); + void pack_y(int,int,int); + void pack_z(int,int,int); + void pack_x_ref(int,int,int); + void pack_y_ref(int,int,int); + void pack_z_ref(int,int,int); +}; + +} // namespace LAMMPS_NS + +#endif + +/* ERROR/WARNING messages: + +E: Cannot find fix store/local + +Fix id cannot be found. + +E: Illegal bond_style command + +Self-explanatory. + +E: Bond style bpm/rotational must include at least one value to output + +Must include at least one bond property to store in fix store/local + +E: Bond style bpm/rotational cannot be used with 3,4-body interactions + +No angle, dihedral, or improper styles can be defined when using +bond style bpm/rotational. + +E: Bond style bpm/rotational cannot be used with atom style template + +This bond style can change the bond topology which is not +allowed with this atom style. + +*/ diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp index 4268655ca1..adc4f596c5 100644 --- a/src/BPM/bond_bpm_rotational.cpp +++ b/src/BPM/bond_bpm_rotational.cpp @@ -1,721 +1,715 @@ -// clang-format off -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "bond_bpm_rotational.h" - -#include "atom.h" -#include "comm.h" -#include "domain.h" -#include "error.h" -#include "fix_bond_history.h" -#include "force.h" -#include "math_const.h" -#include "math_extra.h" -#include "memory.h" -#include "modify.h" -#include "neighbor.h" -#include "pair.h" - -#include -#include -#include - -#define EPSILON 1e-10 - -using namespace LAMMPS_NS; -using namespace MathExtra; - -/* ---------------------------------------------------------------------- */ - -BondBPMRotational::BondBPMRotational(LAMMPS *lmp) : BondBPM(lmp) -{ - partial_flag = 1; - fix_bond_history = nullptr; -} - -/* ---------------------------------------------------------------------- */ - -BondBPMRotational::~BondBPMRotational() -{ - if(fix_bond_history) modify->delete_fix("BOND_HISTORY_BPM_ROTATIONAL"); - - if (allocated) { - memory->destroy(setflag); - memory->destroy(Kr); - memory->destroy(Ks); - memory->destroy(Kt); - memory->destroy(Kb); - memory->destroy(Fcr); - memory->destroy(Fcs); - memory->destroy(Tct); - memory->destroy(Tcb); - memory->destroy(gnorm); - memory->destroy(gslide); - memory->destroy(groll); - memory->destroy(gtwist); - } -} - -/* ---------------------------------------------------------------------- */ - -double BondBPMRotational::acos_limit(double c) -{ - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - return acos(c); -} - -/* ---------------------------------------------------------------------- - Store data for a single bond - if bond added after LAMMPS init (e.g. pour) -------------------------------------------------------------------------- */ - -double BondBPMRotational::store_bond(int n,int i,int j) -{ - int m,k; - double delx, dely, delz, r, rinv; - double **x = atom->x; - tagint *tag = atom->tag; - double **bondstore = fix_bond_history->bondstore; - - if (tag[i] < tag[j]) { - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - } else { - delx = x[j][0] - x[i][0]; - dely = x[j][1] - x[i][1]; - delz = x[j][2] - x[i][2]; - } - - r = sqrt(delx*delx + dely*dely + delz*delz); - rinv = 1.0/r; - - bondstore[n][0] = r; - bondstore[n][1] = delx*rinv; - bondstore[n][2] = dely*rinv; - bondstore[n][3] = delz*rinv; - - if (i < atom->nlocal) { - for (m = 0; m < atom->num_bond[i]; m ++) { - if (atom->bond_atom[i][m] == tag[j]) { - fix_bond_history->update_atom_value(i, m, 0, r); - fix_bond_history->update_atom_value(i, m, 1, delx*rinv); - fix_bond_history->update_atom_value(i, m, 2, dely*rinv); - fix_bond_history->update_atom_value(i, m, 3, delz*rinv); - } - } - } - - if (j < atom->nlocal) { - for (m = 0; m < atom->num_bond[j]; m ++) { - if (atom->bond_atom[j][m] == tag[i]) { - fix_bond_history->update_atom_value(j, m, 0, r); - fix_bond_history->update_atom_value(j, m, 1, delx*rinv); - fix_bond_history->update_atom_value(j, m, 2, dely*rinv); - fix_bond_history->update_atom_value(j, m, 3, delz*rinv); - } - } - } - - return r; -} - -/* ---------------------------------------------------------------------- - Store data for all bonds called once -------------------------------------------------------------------------- */ - -void BondBPMRotational::store_data() -{ - int i, j, m, type; - double delx, dely, delz, r, rinv; - double **x = atom->x; - int **bond_type = atom->bond_type; - tagint *tag = atom->tag; - - for (i = 0; i < atom->nlocal; i ++) { - for (m = 0; m < atom->num_bond[i]; m ++) { - type = bond_type[i][m]; - - //Skip if bond was turned off - if(type < 0) - continue; - - // map to find index n for tag - j = atom->map(atom->bond_atom[i][m]); - if(j == -1) error->one(FLERR, "Atom missing in BPM bond"); - - // Save orientation as pointing towards small tag - if(tag[i] < tag[j]){ - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - } else { - delx = x[j][0] - x[i][0]; - dely = x[j][1] - x[i][1]; - delz = x[j][2] - x[i][2]; - } - - // Get closest image in case bonded with ghost - domain->minimum_image(delx, dely, delz); - - r = sqrt(delx*delx + dely*dely + delz*delz); - rinv = 1.0/r; - fix_bond_history->update_atom_value(i, m, 0, r); - fix_bond_history->update_atom_value(i, m, 1, delx*rinv); - fix_bond_history->update_atom_value(i, m, 2, dely*rinv); - fix_bond_history->update_atom_value(i, m, 3, delz*rinv); - } - } - - fix_bond_history->post_neighbor(); -} - -/* ---------------------------------------------------------------------- */ - -void BondBPMRotational::compute(int eflag, int vflag) -{ - - if (! fix_bond_history->stored_flag) { - fix_bond_history->stored_flag = true; - store_data(); - } - - int i1,i2,itmp,m,n,type,itype,jtype; - double evdwl,fpair,rsq,ebond; - double q1[4], q2[4], r[3], r0[3]; - double r0_mag, r_mag, r_mag_inv, Fr, Fs_mag; - double Tt_mag, Tb_mag; - double force1on2[3], torque1on2[3], torque2on1[3]; - double breaking, smooth; - double rhat[3], wn1[3], wn2[3], wxn1[3], wxn2[3], vroll[3]; - double w1dotr, w2dotr, v1dotr, v2dotr; - double vn1[3], vn2[3], vt1[3], vt2[3], tmp[3], s1[3], s2[3], tdamp[3]; - double tor1, tor2, tor3, fs1, fs2, fs3; - - double q2inv[4], rb[3], rb_x_r0[3], s[3], t[3], Fs[3]; - double q21[4], qp21[4], Tbp[3], Ttp[3]; - double Tsp[3], Fsp[3], Tt[3], Tb[3], Ts[3], F_rot[3], T_rot[3]; - double mq[4], mqinv[4], Ttmp[3], Ftmp[3], qtmp[4]; - double r0_dot_rb, gamma, c, psi, theta, sin_phi, cos_phi, temp; - double mag_in_plane, mag_out_plane; - - ev_init(eflag,vflag); - - if (vflag_global == 2) - force->pair->vflag_either = force->pair->vflag_global = 1; - - double **cutsq = force->pair->cutsq; - double **x = atom->x; - double **v = atom->v; - double **omega = atom->omega; - double **f = atom->f; - double **torque = atom->torque; - double *radius = atom->radius; - double **quat = atom->quat; - tagint *tag = atom->tag; - int **bondlist = neighbor->bondlist; - int nbondlist = neighbor->nbondlist; - int nlocal = atom->nlocal; - int newton_bond = force->newton_bond; - - double **bondstore = fix_bond_history->bondstore; - - for (n = 0; n < nbondlist; n++) { - - // skip bond if already broken - if (bondlist[n][2] <= 0) continue; - - i1 = bondlist[n][0]; - i2 = bondlist[n][1]; - type = bondlist[n][2]; - r0_mag = bondstore[n][0]; - - // Ensure pair is always ordered such that r0 points in - // a consistent direction and to ensure numerical operations - // are identical to minimize the possibility that a bond straddling - // an mpi grid (newton off) doesn't break on one proc but not the other - if (tag[i2] < tag[i1]) { - itmp = i1; - i1 = i2; - i2 = itmp; - } - - // If bond hasn't been set - should be initialized to zero - if (r0_mag < EPSILON || std::isnan(r0_mag)) - r0_mag = store_bond(n,i1,i2); - - r0[0] = bondstore[n][1]; - r0[1] = bondstore[n][2]; - r0[2] = bondstore[n][3]; - MathExtra::scale3(r0_mag, r0); - - q1[0] = quat[i1][0]; - q1[1] = quat[i1][1]; - q1[2] = quat[i1][2]; - q1[3] = quat[i1][3]; - - q2[0] = quat[i2][0]; - q2[1] = quat[i2][1]; - q2[2] = quat[i2][2]; - q2[3] = quat[i2][3]; - - // Note this is the reverse of Mora & Wang - MathExtra::sub3(x[i1], x[i2], r); - - rsq = MathExtra::lensq3(r); - r_mag = sqrt(rsq); - r_mag_inv = 1.0/r_mag; - MathExtra::scale3(r_mag_inv, r, rhat); - - // ------------------------------------------------------// - // Calculate forces using formulation in: - // 1) Y. Wang Acta Geotechnica 2009 - // 2) P. Mora & Y. Wang Advances in Geomcomputing 2009 - // ------------------------------------------------------// - - // Calculate normal forces, rb = bond vector in particle 1's frame - MathExtra::qconjugate(q2, q2inv); - MathExtra::quatrotvec(q2inv, r, rb); - Fr = Kr[type]*(r_mag - r0_mag); - - MathExtra::scale3(Fr*r_mag_inv, rb, F_rot); - - // Calculate forces due to tangential displacements (no rotation) - r0_dot_rb = dot3(r0, rb); - c = r0_dot_rb*r_mag_inv/r0_mag; - gamma = acos_limit(c); - - MathExtra::cross3(rb, r0, rb_x_r0); - MathExtra::cross3(rb, rb_x_r0, s); - MathExtra::norm3(s); - - MathExtra::scale3(Ks[type]*r_mag*gamma, s, Fs); - - // Calculate torque due to tangential displacements - MathExtra::cross3(r0, rb, t); - MathExtra::norm3(t); - - MathExtra::scale3(0.5*r_mag*Ks[type]*r_mag*gamma, t, Ts); - - // Relative rotation force/torque - // Use representation of X'Y'Z' rotations from Wang, Mora 2009 - temp = r_mag + rb[2]; - if (temp < 0.0) temp = 0.0; - mq[0] = sqrt(2)*0.5*sqrt(temp*r_mag_inv); - - temp = sqrt(rb[0]*rb[0]+rb[1]*rb[1]); - if (temp != 0.0) { - mq[1] = -sqrt(2)*0.5/temp; - temp = r_mag - rb[2]; - if (temp < 0.0) temp = 0.0; - mq[1] *= sqrt(temp*r_mag_inv); - mq[2] = -mq[1]; - mq[1] *= rb[1]; - mq[2] *= rb[0]; - } else { - // If aligned along z axis, x,y terms zero (r_mag-rb[2] = 0) - mq[1] = 0.0; - mq[2] = 0.0; - } - mq[3] = 0.0; - - // qp21 = opposite of r^\circ_21 in Wang - // q21 = opposite of r_21 in Wang - MathExtra::quatquat(q2inv, q1, qp21); - MathExtra::qconjugate(mq, mqinv); - MathExtra::quatquat(mqinv,qp21,qtmp); - MathExtra::quatquat(qtmp,mq,q21); - - temp = sqrt(q21[0]*q21[0] + q21[3]*q21[3]); - if (temp != 0.0) { - c = q21[0]/temp; - psi = 2.0*acos_limit(c); - } else { - c = 0.0; - psi = 0.0; - } - - // Map negative rotations - if (q21[3] < 0.0) // sin = q21[3]/temp - psi = -psi; - - if (q21[3] == 0.0) - psi = 0.0; - - c = q21[0]*q21[0] - q21[1]*q21[1] - q21[2]*q21[2] + q21[3]*q21[3]; - theta = acos_limit(c); - - // Separately calculte magnitude of quaternion in x-y and out of x-y planes - // to avoid dividing by zero - mag_out_plane = (q21[0]*q21[0] + q21[3]*q21[3]); - mag_in_plane = (q21[1]*q21[1] + q21[2]*q21[2]); - - if (mag_in_plane == 0.0) { - // No rotation => no bending/shear torque or extra shear force - // achieve by setting cos/sin = 0 - cos_phi = 0.0; - sin_phi = 0.0; - } else if (mag_out_plane == 0.0) { - // Calculate angle in plane - cos_phi = q21[2]/sqrt(mag_in_plane); - sin_phi = -q21[1]/sqrt(mag_in_plane); - } else { - // Default equations in Mora, Wang 2009 - cos_phi = q21[1]*q21[3] + q21[0]*q21[2]; - sin_phi = q21[2]*q21[3] - q21[0]*q21[1]; - - cos_phi /= sqrt(mag_out_plane*mag_in_plane); - sin_phi /= sqrt(mag_out_plane*mag_in_plane); - } - - Tbp[0] = -Kb[type]*theta*sin_phi; - Tbp[1] = Kb[type]*theta*cos_phi; - Tbp[2] = 0.0; - - Ttp[0] = 0.0; - Ttp[1] = 0.0; - Ttp[2] = Kt[type]*psi; - - Fsp[0] = -0.5*Ks[type]*r_mag*theta*cos_phi; - Fsp[1] = -0.5*Ks[type]*r_mag*theta*sin_phi; - Fsp[2] = 0.0; - - Tsp[0] = 0.25*Ks[type]*r_mag*r_mag*theta*sin_phi; - Tsp[1] = -0.25*Ks[type]*r_mag*r_mag*theta*cos_phi; - Tsp[2] = 0.0; - - // Rotate forces/torques back to 1st particle's frame - - MathExtra::quatrotvec(mq, Fsp, Ftmp); - MathExtra::quatrotvec(mq, Tsp, Ttmp); - for (m = 0; m < 3; m++) { - Fs[m] += Ftmp[m]; - Ts[m] += Ttmp[m]; - } - - MathExtra::quatrotvec(mq, Tbp, Tb); - MathExtra::quatrotvec(mq, Ttp, Tt); - - // Sum forces and calculate magnitudes - F_rot[0] += Fs[0]; - F_rot[1] += Fs[1]; - F_rot[2] += Fs[2]; - MathExtra::quatrotvec(q2, F_rot, force1on2); - - T_rot[0] = Ts[0] + Tt[0] + Tb[0]; - T_rot[1] = Ts[1] + Tt[1] + Tb[1]; - T_rot[2] = Ts[2] + Tt[2] + Tb[2]; - MathExtra::quatrotvec(q2, T_rot, torque1on2); - - T_rot[0] = Ts[0] - Tt[0] - Tb[0]; - T_rot[1] = Ts[1] - Tt[1] - Tb[1]; - T_rot[2] = Ts[2] - Tt[2] - Tb[2]; - MathExtra::quatrotvec(q2, T_rot, torque2on1); - - Fs_mag = MathExtra::len3(Fs); - Tt_mag = MathExtra::len3(Tt); - Tb_mag = MathExtra::len3(Tb); - - // ------------------------------------------------------// - // Check if bond breaks - // ------------------------------------------------------// - - breaking = Fs_mag/Fcs[type] + Tb_mag/Tcb[type] + Tt_mag/Tct[type]; - if (Fr > 0.0) breaking += Fr/Fcr[type]; - - if (breaking >= 1.0) { - bondlist[n][2] = 0; - process_broken(i1, i2); - continue; - } - - smooth = breaking*breaking; - smooth = 1.0 - smooth*smooth; - - // Not actual energy, just a handy metric - if (eflag) ebond = -smooth; - - // ------------------------------------------------------// - // Calculate damping using formulation in - // Y. Wang, F. Alonso-Marroquin, W. Guo 2015 - // ------------------------------------------------------// - // Note: n points towards 1 vs pointing towards 2 - - // Damp normal velocity difference - v1dotr = MathExtra::dot3(v[i1],rhat); - v2dotr = MathExtra::dot3(v[i2],rhat); - - MathExtra::scale3(v1dotr, rhat, vn1); - MathExtra::scale3(v2dotr, rhat, vn2); - - MathExtra::sub3(vn1, vn2, tmp); - MathExtra::scale3(gnorm[type], tmp); - MathExtra::add3(force1on2, tmp, force1on2); - - // Damp tangential objective velocities - MathExtra::sub3(v[i1], vn1, vt1); - MathExtra::sub3(v[i2], vn2, vt2); - - MathExtra::sub3(vt2, vt1, tmp); - MathExtra::scale3(-0.5, tmp); - - MathExtra::cross3(omega[i1], r, s1); - MathExtra::scale3(0.5, s1); - MathExtra::sub3(s1, tmp, s1); // Eq 12 - - MathExtra::cross3(omega[i2], r, s2); - MathExtra::scale3(-0.5,s2); - MathExtra::add3(s2, tmp, s2); // Eq 13 - MathExtra::scale3(-0.5,s2); - - MathExtra::sub3(s1, s2, tmp); - MathExtra::scale3(gslide[type], tmp); - MathExtra::add3(force1on2, tmp, force1on2); - - // Apply corresponding torque - MathExtra::cross3(r,tmp,tdamp); - MathExtra::scale3(-0.5, tdamp); // 0.5*r points from particle 2 to midpoint - MathExtra::add3(torque1on2, tdamp, torque1on2); - MathExtra::add3(torque2on1, tdamp, torque2on1); - - // Damp rolling - MathExtra::cross3(omega[i1], rhat, wxn1); - MathExtra::cross3(omega[i2], rhat, wxn2); - MathExtra::sub3(wxn1, wxn2, vroll); // Eq. 31 - MathExtra::cross3(r, vroll, tdamp); - - MathExtra::scale3(0.5*groll[type], tdamp); - MathExtra::add3(torque1on2, tdamp, torque1on2); - MathExtra::scale3(-1.0, tdamp); - MathExtra::add3(torque2on1, tdamp, torque2on1); - - // Damp twist - w1dotr = MathExtra::dot3(omega[i1],rhat); - w2dotr = MathExtra::dot3(omega[i2],rhat); - - MathExtra::scale3(w1dotr, rhat, wn1); - MathExtra::scale3(w2dotr, rhat, wn2); - - MathExtra::sub3(wn1, wn2, tdamp); // Eq. 38 - MathExtra::scale3(0.5*gtwist[type], tdamp); - MathExtra::add3(torque1on2, tdamp, torque1on2); - MathExtra::scale3(-1.0, tdamp); - MathExtra::add3(torque2on1, tdamp, torque2on1); - - // ------------------------------------------------------// - // Apply forces and torques to particles - // ------------------------------------------------------// - - if (newton_bond || i1 < nlocal) { - f[i1][0] -= force1on2[0]*smooth; - f[i1][1] -= force1on2[1]*smooth; - f[i1][2] -= force1on2[2]*smooth; - - torque[i1][0] += torque2on1[0]*smooth; - torque[i1][1] += torque2on1[1]*smooth; - torque[i1][2] += torque2on1[2]*smooth; - } - - if (newton_bond || i2 < nlocal) { - f[i2][0] += force1on2[0]*smooth; - f[i2][1] += force1on2[1]*smooth; - f[i2][2] += force1on2[2]*smooth; - - torque[i2][0] += torque1on2[0]*smooth; - torque[i2][1] += torque1on2[1]*smooth; - torque[i2][2] += torque1on2[2]*smooth; - } - - if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,Fr*smooth,r[0],r[1],r[2]); - } -} - -/* ---------------------------------------------------------------------- */ - -void BondBPMRotational::allocate() -{ - allocated = 1; - int n = atom->nbondtypes; - - memory->create(Kr,n+1,"bond:Kr"); - memory->create(Ks,n+1,"bond:Ks"); - memory->create(Kt,n+1,"bond:Kt"); - memory->create(Kb,n+1,"bond:Kb"); - memory->create(Fcr,n+1,"bond:Fcr"); - memory->create(Fcs,n+1,"bond:Fcs"); - memory->create(Tct,n+1,"bond:Tct"); - memory->create(Tcb,n+1,"bond:Tcb"); - memory->create(gnorm,n+1,"bond:gnorm"); - memory->create(gslide,n+1,"bond:gslide"); - memory->create(groll,n+1,"bond:groll"); - memory->create(gtwist,n+1,"bond:gtwist"); - - memory->create(setflag,n+1,"bond:setflag"); - for (int i = 1; i <= n; i++) setflag[i] = 0; -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more types -------------------------------------------------------------------------- */ - -void BondBPMRotational::coeff(int narg, char **arg) -{ - if (narg != 13) error->all(FLERR,"Incorrect args for bond coefficients"); - if (!allocated) allocate(); - - int ilo,ihi; - utils::bounds(FLERR,arg[0],1,atom->nbondtypes,ilo,ihi,error); - - double Kr_one = utils::numeric(FLERR,arg[1],false,lmp); - double Ks_one = utils::numeric(FLERR,arg[2],false,lmp); - double Kt_one = utils::numeric(FLERR,arg[3],false,lmp); - double Kb_one = utils::numeric(FLERR,arg[4],false,lmp); - double Fcr_one = utils::numeric(FLERR,arg[5],false,lmp); - double Fcs_one = utils::numeric(FLERR,arg[6],false,lmp); - double Tct_one = utils::numeric(FLERR,arg[7],false,lmp); - double Tcb_one = utils::numeric(FLERR,arg[8],false,lmp); - double gnorm_one = utils::numeric(FLERR,arg[9],false,lmp); - double gslide_one = utils::numeric(FLERR,arg[10],false,lmp); - double groll_one = utils::numeric(FLERR,arg[11],false,lmp); - double gtwist_one = utils::numeric(FLERR,arg[12],false,lmp); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - Kr[i] = Kr_one; - Ks[i] = Ks_one; - Kt[i] = Kt_one; - Kb[i] = Kb_one; - Fcr[i] = Fcr_one; - Fcs[i] = Fcs_one; - Tct[i] = Tct_one; - Tcb[i] = Tcb_one; - gnorm[i] = gnorm_one; - gslide[i] = gslide_one; - groll[i] = groll_one; - gtwist[i] = gtwist_one; - setflag[i] = 1; - count++; - - if (Fcr[i]/Kr[i] > max_stretch) max_stretch = Fcr[i]/Kr[i]; - } - - if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients"); -} - -/* ---------------------------------------------------------------------- - check if pair defined and special_bond settings are valid -------------------------------------------------------------------------- */ - -void BondBPMRotational::init_style() -{ - BondBPM::init_style(); - - if (!atom->quat_flag || !atom->sphere_flag) - error->all(FLERR,"Bond bpm/rotational requires atom style sphere/bpm"); - if (comm->ghost_velocity == 0) - error->all(FLERR,"Bond bpm/rotational requires ghost atoms store velocity"); - - if(domain->dimension == 2) - error->warning(FLERR, "Bond style bpm/rotational not intended for 2d use"); - - if (!fix_bond_history) - fix_bond_history = (FixBondHistory *) modify->add_fix( - "BOND_HISTORY_BPM_ROTATIONAL all BOND_HISTORY 0 4"); -} - -/* ---------------------------------------------------------------------- - proc 0 writes out coeffs to restart file -------------------------------------------------------------------------- */ - -void BondBPMRotational::write_restart(FILE *fp) -{ - fwrite(&Kr[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&Ks[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&Kt[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&Kb[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&Fcr[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&Fcs[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&Tct[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&Tcb[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&gnorm[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&gslide[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&groll[1],sizeof(double),atom->nbondtypes,fp); - fwrite(>wist[1],sizeof(double),atom->nbondtypes,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads coeffs from restart file, bcasts them -------------------------------------------------------------------------- */ - -void BondBPMRotational::read_restart(FILE *fp) -{ - allocate(); - - if (comm->me == 0) { - utils::sfread(FLERR,&Kr[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,&Ks[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,&Kt[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,&Kb[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,&Fcr[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,&Fcs[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,&Tct[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,&Tcb[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,&gnorm[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,&gslide[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,&groll[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - utils::sfread(FLERR,>wist[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); - } - MPI_Bcast(&Kr[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&Ks[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&Kt[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&Kb[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&Fcr[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&Fcs[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&Tct[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&Tcb[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&gnorm[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&gslide[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&groll[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(>wist[1],atom->nbondtypes,MPI_DOUBLE,0,world); - - for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to data file -------------------------------------------------------------------------- */ - -void BondBPMRotational::write_data(FILE *fp) -{ - for (int i = 1; i <= atom->nbondtypes; i++) - fprintf(fp,"%d %g %g %g %g %g %g %g %g %g %g %g %g\n", - i,Kr[i],Ks[i],Kt[i],Kb[i],Fcr[i], Fcs[i], Tct[i], - Tcb[i], gnorm[i], gslide[i], groll[i], gtwist[i]); -} - -/* ---------------------------------------------------------------------- */ - -double BondBPMRotational::single(int type, double rsq, int i, int j, - double &fforce) -{ - // Not yet enabled - if (type <= 0) return 0.0; - - //double r0; - //for (int n = 0; n < atom->num_bond[i]; n ++) { - // if (atom->bond_atom[i][n] == atom->tag[j]) { - // r0 = fix_bond_history->get_atom_value(i, n, 0); - // } - //} -} +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "bond_bpm_rotational.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "fix_bond_history.h" +#include "force.h" +#include "math_const.h" +#include "math_extra.h" +#include "memory.h" +#include "modify.h" +#include "neighbor.h" +#include "pair.h" + +#include +#include +#include + +#define EPSILON 1e-10 + +using namespace LAMMPS_NS; +using namespace MathExtra; + +/* ---------------------------------------------------------------------- */ + +BondBPMRotational::BondBPMRotational(LAMMPS *lmp) : BondBPM(lmp) +{ + partial_flag = 1; +} + +/* ---------------------------------------------------------------------- */ + +BondBPMRotational::~BondBPMRotational() +{ + if (fix_bond_history) modify->delete_fix("BOND_HISTORY_BPM_ROTATIONAL"); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(Kr); + memory->destroy(Ks); + memory->destroy(Kt); + memory->destroy(Kb); + memory->destroy(Fcr); + memory->destroy(Fcs); + memory->destroy(Tct); + memory->destroy(Tcb); + memory->destroy(gnorm); + memory->destroy(gslide); + memory->destroy(groll); + memory->destroy(gtwist); + } +} + +/* ---------------------------------------------------------------------- */ + +double BondBPMRotational::acos_limit(double c) +{ + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + return acos(c); +} + +/* ---------------------------------------------------------------------- + Store data for a single bond - if bond added after LAMMPS init (e.g. pour) +------------------------------------------------------------------------- */ + +double BondBPMRotational::store_bond(int n,int i,int j) +{ + int m,k; + double delx, dely, delz, r, rinv; + double **x = atom->x; + tagint *tag = atom->tag; + double **bondstore = fix_bond_history->bondstore; + + if (tag[i] < tag[j]) { + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + } else { + delx = x[j][0] - x[i][0]; + dely = x[j][1] - x[i][1]; + delz = x[j][2] - x[i][2]; + } + + r = sqrt(delx*delx + dely*dely + delz*delz); + rinv = 1.0/r; + + bondstore[n][0] = r; + bondstore[n][1] = delx*rinv; + bondstore[n][2] = dely*rinv; + bondstore[n][3] = delz*rinv; + + if (i < atom->nlocal) { + for (m = 0; m < atom->num_bond[i]; m ++) { + if (atom->bond_atom[i][m] == tag[j]) { + fix_bond_history->update_atom_value(i, m, 0, r); + fix_bond_history->update_atom_value(i, m, 1, delx*rinv); + fix_bond_history->update_atom_value(i, m, 2, dely*rinv); + fix_bond_history->update_atom_value(i, m, 3, delz*rinv); + } + } + } + + if (j < atom->nlocal) { + for (m = 0; m < atom->num_bond[j]; m ++) { + if (atom->bond_atom[j][m] == tag[i]) { + fix_bond_history->update_atom_value(j, m, 0, r); + fix_bond_history->update_atom_value(j, m, 1, delx*rinv); + fix_bond_history->update_atom_value(j, m, 2, dely*rinv); + fix_bond_history->update_atom_value(j, m, 3, delz*rinv); + } + } + } + + return r; +} + +/* ---------------------------------------------------------------------- + Store data for all bonds called once +------------------------------------------------------------------------- */ + +void BondBPMRotational::store_data() +{ + int i, j, m, type; + double delx, dely, delz, r, rinv; + double **x = atom->x; + int **bond_type = atom->bond_type; + tagint *tag = atom->tag; + + for (i = 0; i < atom->nlocal; i ++) { + for (m = 0; m < atom->num_bond[i]; m ++) { + type = bond_type[i][m]; + + //Skip if bond was turned off + if(type < 0) + continue; + + // map to find index n for tag + j = atom->map(atom->bond_atom[i][m]); + if(j == -1) error->one(FLERR, "Atom missing in BPM bond"); + + // Save orientation as pointing towards small tag + if(tag[i] < tag[j]){ + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + } else { + delx = x[j][0] - x[i][0]; + dely = x[j][1] - x[i][1]; + delz = x[j][2] - x[i][2]; + } + + // Get closest image in case bonded with ghost + domain->minimum_image(delx, dely, delz); + + r = sqrt(delx*delx + dely*dely + delz*delz); + rinv = 1.0/r; + fix_bond_history->update_atom_value(i, m, 0, r); + fix_bond_history->update_atom_value(i, m, 1, delx*rinv); + fix_bond_history->update_atom_value(i, m, 2, dely*rinv); + fix_bond_history->update_atom_value(i, m, 3, delz*rinv); + } + } + + fix_bond_history->post_neighbor(); +} + +/* ---------------------------------------------------------------------- */ + +void BondBPMRotational::compute(int eflag, int vflag) +{ + + if (! fix_bond_history->stored_flag) { + fix_bond_history->stored_flag = true; + store_data(); + } + + int i1,i2,itmp,m,n,type,itype,jtype; + double q1[4], q2[4], r[3], r0[3]; + double rsq, r0_mag, r_mag, r_mag_inv, Fr, Fs_mag; + double Tt_mag, Tb_mag, breaking, smooth; + double force1on2[3], torque1on2[3], torque2on1[3]; + double rhat[3], wn1[3], wn2[3], wxn1[3], wxn2[3], vroll[3]; + double w1dotr, w2dotr, v1dotr, v2dotr; + double vn1[3], vn2[3], vt1[3], vt2[3], tmp[3], s1[3], s2[3], tdamp[3]; + double tor1, tor2, tor3, fs1, fs2, fs3; + + double q2inv[4], rb[3], rb_x_r0[3], s[3], t[3], Fs[3]; + double q21[4], qp21[4], Tbp[3], Ttp[3]; + double Tsp[3], Fsp[3], Tt[3], Tb[3], Ts[3], F_rot[3], T_rot[3]; + double mq[4], mqinv[4], Ttmp[3], Ftmp[3], qtmp[4]; + double r0_dot_rb, gamma, c, psi, theta, sin_phi, cos_phi, temp; + double mag_in_plane, mag_out_plane; + + ev_init(eflag,vflag); + + if (vflag_global == 2) + force->pair->vflag_either = force->pair->vflag_global = 1; + + double **cutsq = force->pair->cutsq; + double **x = atom->x; + double **v = atom->v; + double **omega = atom->omega; + double **f = atom->f; + double **torque = atom->torque; + double *radius = atom->radius; + double **quat = atom->quat; + tagint *tag = atom->tag; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + double **bondstore = fix_bond_history->bondstore; + + for (n = 0; n < nbondlist; n++) { + + // skip bond if already broken + if (bondlist[n][2] <= 0) continue; + + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + type = bondlist[n][2]; + r0_mag = bondstore[n][0]; + + // Ensure pair is always ordered such that r0 points in + // a consistent direction and to ensure numerical operations + // are identical to minimize the possibility that a bond straddling + // an mpi grid (newton off) doesn't break on one proc but not the other + if (tag[i2] < tag[i1]) { + itmp = i1; + i1 = i2; + i2 = itmp; + } + + // If bond hasn't been set - should be initialized to zero + if (r0_mag < EPSILON || std::isnan(r0_mag)) + r0_mag = store_bond(n,i1,i2); + + r0[0] = bondstore[n][1]; + r0[1] = bondstore[n][2]; + r0[2] = bondstore[n][3]; + MathExtra::scale3(r0_mag, r0); + + q1[0] = quat[i1][0]; + q1[1] = quat[i1][1]; + q1[2] = quat[i1][2]; + q1[3] = quat[i1][3]; + + q2[0] = quat[i2][0]; + q2[1] = quat[i2][1]; + q2[2] = quat[i2][2]; + q2[3] = quat[i2][3]; + + // Note this is the reverse of Mora & Wang + MathExtra::sub3(x[i1], x[i2], r); + + rsq = MathExtra::lensq3(r); + r_mag = sqrt(rsq); + r_mag_inv = 1.0/r_mag; + MathExtra::scale3(r_mag_inv, r, rhat); + + // ------------------------------------------------------// + // Calculate forces using formulation in: + // 1) Y. Wang Acta Geotechnica 2009 + // 2) P. Mora & Y. Wang Advances in Geomcomputing 2009 + // ------------------------------------------------------// + + // Calculate normal forces, rb = bond vector in particle 1's frame + MathExtra::qconjugate(q2, q2inv); + MathExtra::quatrotvec(q2inv, r, rb); + Fr = Kr[type]*(r_mag - r0_mag); + + MathExtra::scale3(Fr*r_mag_inv, rb, F_rot); + + // Calculate forces due to tangential displacements (no rotation) + r0_dot_rb = dot3(r0, rb); + c = r0_dot_rb*r_mag_inv/r0_mag; + gamma = acos_limit(c); + + MathExtra::cross3(rb, r0, rb_x_r0); + MathExtra::cross3(rb, rb_x_r0, s); + MathExtra::norm3(s); + + MathExtra::scale3(Ks[type]*r_mag*gamma, s, Fs); + + // Calculate torque due to tangential displacements + MathExtra::cross3(r0, rb, t); + MathExtra::norm3(t); + + MathExtra::scale3(0.5*r_mag*Ks[type]*r_mag*gamma, t, Ts); + + // Relative rotation force/torque + // Use representation of X'Y'Z' rotations from Wang, Mora 2009 + temp = r_mag + rb[2]; + if (temp < 0.0) temp = 0.0; + mq[0] = sqrt(2)*0.5*sqrt(temp*r_mag_inv); + + temp = sqrt(rb[0]*rb[0]+rb[1]*rb[1]); + if (temp != 0.0) { + mq[1] = -sqrt(2)*0.5/temp; + temp = r_mag - rb[2]; + if (temp < 0.0) temp = 0.0; + mq[1] *= sqrt(temp*r_mag_inv); + mq[2] = -mq[1]; + mq[1] *= rb[1]; + mq[2] *= rb[0]; + } else { + // If aligned along z axis, x,y terms zero (r_mag-rb[2] = 0) + mq[1] = 0.0; + mq[2] = 0.0; + } + mq[3] = 0.0; + + // qp21 = opposite of r^\circ_21 in Wang + // q21 = opposite of r_21 in Wang + MathExtra::quatquat(q2inv, q1, qp21); + MathExtra::qconjugate(mq, mqinv); + MathExtra::quatquat(mqinv,qp21,qtmp); + MathExtra::quatquat(qtmp,mq,q21); + + temp = sqrt(q21[0]*q21[0] + q21[3]*q21[3]); + if (temp != 0.0) { + c = q21[0]/temp; + psi = 2.0*acos_limit(c); + } else { + c = 0.0; + psi = 0.0; + } + + // Map negative rotations + if (q21[3] < 0.0) // sin = q21[3]/temp + psi = -psi; + + if (q21[3] == 0.0) + psi = 0.0; + + c = q21[0]*q21[0] - q21[1]*q21[1] - q21[2]*q21[2] + q21[3]*q21[3]; + theta = acos_limit(c); + + // Separately calculte magnitude of quaternion in x-y and out of x-y planes + // to avoid dividing by zero + mag_out_plane = (q21[0]*q21[0] + q21[3]*q21[3]); + mag_in_plane = (q21[1]*q21[1] + q21[2]*q21[2]); + + if (mag_in_plane == 0.0) { + // No rotation => no bending/shear torque or extra shear force + // achieve by setting cos/sin = 0 + cos_phi = 0.0; + sin_phi = 0.0; + } else if (mag_out_plane == 0.0) { + // Calculate angle in plane + cos_phi = q21[2]/sqrt(mag_in_plane); + sin_phi = -q21[1]/sqrt(mag_in_plane); + } else { + // Default equations in Mora, Wang 2009 + cos_phi = q21[1]*q21[3] + q21[0]*q21[2]; + sin_phi = q21[2]*q21[3] - q21[0]*q21[1]; + + cos_phi /= sqrt(mag_out_plane*mag_in_plane); + sin_phi /= sqrt(mag_out_plane*mag_in_plane); + } + + Tbp[0] = -Kb[type]*theta*sin_phi; + Tbp[1] = Kb[type]*theta*cos_phi; + Tbp[2] = 0.0; + + Ttp[0] = 0.0; + Ttp[1] = 0.0; + Ttp[2] = Kt[type]*psi; + + Fsp[0] = -0.5*Ks[type]*r_mag*theta*cos_phi; + Fsp[1] = -0.5*Ks[type]*r_mag*theta*sin_phi; + Fsp[2] = 0.0; + + Tsp[0] = 0.25*Ks[type]*r_mag*r_mag*theta*sin_phi; + Tsp[1] = -0.25*Ks[type]*r_mag*r_mag*theta*cos_phi; + Tsp[2] = 0.0; + + // Rotate forces/torques back to 1st particle's frame + + MathExtra::quatrotvec(mq, Fsp, Ftmp); + MathExtra::quatrotvec(mq, Tsp, Ttmp); + for (m = 0; m < 3; m++) { + Fs[m] += Ftmp[m]; + Ts[m] += Ttmp[m]; + } + + MathExtra::quatrotvec(mq, Tbp, Tb); + MathExtra::quatrotvec(mq, Ttp, Tt); + + // Sum forces and calculate magnitudes + F_rot[0] += Fs[0]; + F_rot[1] += Fs[1]; + F_rot[2] += Fs[2]; + MathExtra::quatrotvec(q2, F_rot, force1on2); + + T_rot[0] = Ts[0] + Tt[0] + Tb[0]; + T_rot[1] = Ts[1] + Tt[1] + Tb[1]; + T_rot[2] = Ts[2] + Tt[2] + Tb[2]; + MathExtra::quatrotvec(q2, T_rot, torque1on2); + + T_rot[0] = Ts[0] - Tt[0] - Tb[0]; + T_rot[1] = Ts[1] - Tt[1] - Tb[1]; + T_rot[2] = Ts[2] - Tt[2] - Tb[2]; + MathExtra::quatrotvec(q2, T_rot, torque2on1); + + Fs_mag = MathExtra::len3(Fs); + Tt_mag = MathExtra::len3(Tt); + Tb_mag = MathExtra::len3(Tb); + + // ------------------------------------------------------// + // Check if bond breaks + // ------------------------------------------------------// + + breaking = Fr/Fcr[type] + Fs_mag/Fcs[type] + Tb_mag/Tcb[type] + Tt_mag/Tct[type]; + if (breaking < 0.0) breaking = 0.0; + + if (breaking >= 1.0) { + bondlist[n][2] = 0; + process_broken(i1, i2); + continue; + } + + smooth = breaking*breaking; + smooth = 1.0 - smooth*smooth; + + // ------------------------------------------------------// + // Calculate damping using formulation in + // Y. Wang, F. Alonso-Marroquin, W. Guo 2015 + // ------------------------------------------------------// + // Note: n points towards 1 vs pointing towards 2 + + // Damp normal velocity difference + v1dotr = MathExtra::dot3(v[i1],rhat); + v2dotr = MathExtra::dot3(v[i2],rhat); + + MathExtra::scale3(v1dotr, rhat, vn1); + MathExtra::scale3(v2dotr, rhat, vn2); + + MathExtra::sub3(vn1, vn2, tmp); + MathExtra::scale3(gnorm[type], tmp); + MathExtra::add3(force1on2, tmp, force1on2); + + // Damp tangential objective velocities + MathExtra::sub3(v[i1], vn1, vt1); + MathExtra::sub3(v[i2], vn2, vt2); + + MathExtra::sub3(vt2, vt1, tmp); + MathExtra::scale3(-0.5, tmp); + + MathExtra::cross3(omega[i1], r, s1); + MathExtra::scale3(0.5, s1); + MathExtra::sub3(s1, tmp, s1); // Eq 12 + + MathExtra::cross3(omega[i2], r, s2); + MathExtra::scale3(-0.5,s2); + MathExtra::add3(s2, tmp, s2); // Eq 13 + MathExtra::scale3(-0.5,s2); + + MathExtra::sub3(s1, s2, tmp); + MathExtra::scale3(gslide[type], tmp); + MathExtra::add3(force1on2, tmp, force1on2); + + // Apply corresponding torque + MathExtra::cross3(r,tmp,tdamp); + MathExtra::scale3(-0.5, tdamp); // 0.5*r points from particle 2 to midpoint + MathExtra::add3(torque1on2, tdamp, torque1on2); + MathExtra::add3(torque2on1, tdamp, torque2on1); + + // Damp rolling + MathExtra::cross3(omega[i1], rhat, wxn1); + MathExtra::cross3(omega[i2], rhat, wxn2); + MathExtra::sub3(wxn1, wxn2, vroll); // Eq. 31 + MathExtra::cross3(r, vroll, tdamp); + + MathExtra::scale3(0.5*groll[type], tdamp); + MathExtra::add3(torque1on2, tdamp, torque1on2); + MathExtra::scale3(-1.0, tdamp); + MathExtra::add3(torque2on1, tdamp, torque2on1); + + // Damp twist + w1dotr = MathExtra::dot3(omega[i1],rhat); + w2dotr = MathExtra::dot3(omega[i2],rhat); + + MathExtra::scale3(w1dotr, rhat, wn1); + MathExtra::scale3(w2dotr, rhat, wn2); + + MathExtra::sub3(wn1, wn2, tdamp); // Eq. 38 + MathExtra::scale3(0.5*gtwist[type], tdamp); + MathExtra::add3(torque1on2, tdamp, torque1on2); + MathExtra::scale3(-1.0, tdamp); + MathExtra::add3(torque2on1, tdamp, torque2on1); + + // ------------------------------------------------------// + // Apply forces and torques to particles + // ------------------------------------------------------// + + if (newton_bond || i1 < nlocal) { + f[i1][0] -= force1on2[0]*smooth; + f[i1][1] -= force1on2[1]*smooth; + f[i1][2] -= force1on2[2]*smooth; + + torque[i1][0] += torque2on1[0]*smooth; + torque[i1][1] += torque2on1[1]*smooth; + torque[i1][2] += torque2on1[2]*smooth; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] += force1on2[0]*smooth; + f[i2][1] += force1on2[1]*smooth; + f[i2][2] += force1on2[2]*smooth; + + torque[i2][0] += torque1on2[0]*smooth; + torque[i2][1] += torque1on2[1]*smooth; + torque[i2][2] += torque1on2[2]*smooth; + } + + if (evflag) ev_tally(i1,i2,nlocal,newton_bond,0.0,Fr*smooth,r[0],r[1],r[2]); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondBPMRotational::allocate() +{ + allocated = 1; + int n = atom->nbondtypes; + + memory->create(Kr,n+1,"bond:Kr"); + memory->create(Ks,n+1,"bond:Ks"); + memory->create(Kt,n+1,"bond:Kt"); + memory->create(Kb,n+1,"bond:Kb"); + memory->create(Fcr,n+1,"bond:Fcr"); + memory->create(Fcs,n+1,"bond:Fcs"); + memory->create(Tct,n+1,"bond:Tct"); + memory->create(Tcb,n+1,"bond:Tcb"); + memory->create(gnorm,n+1,"bond:gnorm"); + memory->create(gslide,n+1,"bond:gslide"); + memory->create(groll,n+1,"bond:groll"); + memory->create(gtwist,n+1,"bond:gtwist"); + + memory->create(setflag,n+1,"bond:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void BondBPMRotational::coeff(int narg, char **arg) +{ + if (narg != 13) error->all(FLERR,"Incorrect args for bond coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + utils::bounds(FLERR,arg[0],1,atom->nbondtypes,ilo,ihi,error); + + double Kr_one = utils::numeric(FLERR,arg[1],false,lmp); + double Ks_one = utils::numeric(FLERR,arg[2],false,lmp); + double Kt_one = utils::numeric(FLERR,arg[3],false,lmp); + double Kb_one = utils::numeric(FLERR,arg[4],false,lmp); + double Fcr_one = utils::numeric(FLERR,arg[5],false,lmp); + double Fcs_one = utils::numeric(FLERR,arg[6],false,lmp); + double Tct_one = utils::numeric(FLERR,arg[7],false,lmp); + double Tcb_one = utils::numeric(FLERR,arg[8],false,lmp); + double gnorm_one = utils::numeric(FLERR,arg[9],false,lmp); + double gslide_one = utils::numeric(FLERR,arg[10],false,lmp); + double groll_one = utils::numeric(FLERR,arg[11],false,lmp); + double gtwist_one = utils::numeric(FLERR,arg[12],false,lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + Kr[i] = Kr_one; + Ks[i] = Ks_one; + Kt[i] = Kt_one; + Kb[i] = Kb_one; + Fcr[i] = Fcr_one; + Fcs[i] = Fcs_one; + Tct[i] = Tct_one; + Tcb[i] = Tcb_one; + gnorm[i] = gnorm_one; + gslide[i] = gslide_one; + groll[i] = groll_one; + gtwist[i] = gtwist_one; + setflag[i] = 1; + count++; + + if (Fcr[i]/Kr[i] > max_stretch) max_stretch = Fcr[i]/Kr[i]; + } + + if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients"); +} + +/* ---------------------------------------------------------------------- + check if pair defined and special_bond settings are valid +------------------------------------------------------------------------- */ + +void BondBPMRotational::init_style() +{ + BondBPM::init_style(); + + if (!atom->quat_flag || !atom->sphere_flag) + error->all(FLERR,"Bond bpm/rotational requires atom style sphere/bpm"); + if (comm->ghost_velocity == 0) + error->all(FLERR,"Bond bpm/rotational requires ghost atoms store velocity"); + + if(domain->dimension == 2) + error->warning(FLERR, "Bond style bpm/rotational not intended for 2d use"); + + if (!fix_bond_history) + fix_bond_history = (FixBondHistory *) modify->add_fix( + "BOND_HISTORY_BPM_ROTATIONAL all BOND_HISTORY 0 4"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void BondBPMRotational::write_restart(FILE *fp) +{ + fwrite(&Kr[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Ks[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Kt[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Kb[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Fcr[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Fcs[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Tct[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Tcb[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&gnorm[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&gslide[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&groll[1],sizeof(double),atom->nbondtypes,fp); + fwrite(>wist[1],sizeof(double),atom->nbondtypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void BondBPMRotational::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR,&Kr[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Ks[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Kt[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Kb[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Fcr[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Fcs[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Tct[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Tcb[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&gnorm[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&gslide[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&groll[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,>wist[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + } + MPI_Bcast(&Kr[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Ks[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Kt[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Kb[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Fcr[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Fcs[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Tct[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Tcb[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&gnorm[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&gslide[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&groll[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(>wist[1],atom->nbondtypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void BondBPMRotational::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nbondtypes; i++) + fprintf(fp,"%d %g %g %g %g %g %g %g %g %g %g %g %g\n", + i,Kr[i],Ks[i],Kt[i],Kb[i],Fcr[i], Fcs[i], Tct[i], + Tcb[i], gnorm[i], gslide[i], groll[i], gtwist[i]); +} + +/* ---------------------------------------------------------------------- */ + +double BondBPMRotational::single(int type, double rsq, int i, int j, + double &fforce) +{ + // Not yet enabled + if (type <= 0) return 0.0; + + //double r0; + //for (int n = 0; n < atom->num_bond[i]; n ++) { + // if (atom->bond_atom[i][n] == atom->tag[j]) { + // r0 = fix_bond_history->get_atom_value(i, n, 0); + // } + //} +} diff --git a/src/BPM/bond_bpm_rotational.h b/src/BPM/bond_bpm_rotational.h index 838226f0ff..14105f11a0 100644 --- a/src/BPM/bond_bpm_rotational.h +++ b/src/BPM/bond_bpm_rotational.h @@ -44,7 +44,6 @@ class BondBPMRotational : public BondBPM { void allocate(); void store_data(); double store_bond(int, int, int); - class FixBondHistory *fix_bond_history; }; } // namespace LAMMPS_NS diff --git a/src/BPM/fix_bond_history.cpp b/src/BPM/fix_bond_history.cpp index b0f046f43f..d66dec5643 100644 --- a/src/BPM/fix_bond_history.cpp +++ b/src/BPM/fix_bond_history.cpp @@ -291,4 +291,16 @@ void FixBondHistory::set_arrays(int i) for (int m = 0; m < nbond; m++) for (int idata = 0; idata < ndata; idata++) stored[i][m*ndata+idata] = 0.0; +} + +/* ---------------------------------------------------------------------- + Remove all data for row by compressing data +------------------------------------------------------------------------- */ + +void FixBondHistory::delete_bond(int i, int m) +{ + double **stored = atom->darray[index]; + int n = atom->num_bond[i]; + for (int idata = 0; idata < ndata; idata ++) + stored[i][m*ndata+idata] = stored[i][(n-1)*ndata+idata]; } \ No newline at end of file diff --git a/src/BPM/fix_bond_history.h b/src/BPM/fix_bond_history.h index 0219e5de02..63102f123c 100644 --- a/src/BPM/fix_bond_history.h +++ b/src/BPM/fix_bond_history.h @@ -41,6 +41,7 @@ class FixBondHistory : public Fix { void update_atom_value(int, int, int, double); double get_atom_value(int, int, int); + void delete_bond(int,int); double **bondstore; int stored_flag;