diff --git a/examples/USER/gle/data.tip4pf b/examples/USER/gle/data.tip4pf new file mode 100644 index 0000000000..a02ff22c50 --- /dev/null +++ b/examples/USER/gle/data.tip4pf @@ -0,0 +1,412 @@ +LAMMPS Description + + 192 atoms + 128 bonds + 64 angles + + 2 atom types + 1 bond types + 1 angle types + + 0 12.4244 xlo xhi + 0 12.4244 ylo yhi + 0 12.4244 zlo zhi + +Masses + + 1 15.9994 + 2 1.0080 + +Atoms + + 1 1 1 -1.1128 2.7475 12.8896 2.6639 + 2 1 2 0.5564 3.1049 12.7379 1.8224 + 3 1 2 0.5564 1.7943 13.1205 2.5592 + 4 2 1 -1.1128 -0.7565 -1.1504 -1.5704 + 5 2 2 0.5564 -1.0109 -0.3411 -2.0812 + 6 2 2 0.5564 -0.7469 -1.7213 -2.2771 + 7 3 1 -1.1128 5.7663 -0.4112 0.4047 + 8 3 2 0.5564 6.4944 -0.7495 -0.3433 + 9 3 2 0.5564 4.9901 -0.2595 -0.2404 + 10 4 1 -1.1128 12.5178 13.4586 2.5833 + 11 4 2 0.5564 11.9595 12.7227 2.5625 + 12 4 2 0.5564 12.3870 13.9685 1.6020 + 13 5 1 -1.1128 11.5658 1.2126 -3.4671 + 14 5 2 0.5564 10.7292 1.5210 -4.0420 + 15 5 2 0.5564 12.2784 0.9456 -4.0624 + 16 6 1 -1.1128 14.6971 6.1882 -7.7474 + 17 6 2 0.5564 14.4784 5.3902 -7.2900 + 18 6 2 0.5564 15.2258 6.8450 -7.3017 + 19 7 1 -1.1128 -2.6192 10.3928 -0.0001 + 20 7 2 0.5564 -1.8647 10.9848 -0.5571 + 21 7 2 0.5564 -2.7265 9.2890 -0.0662 + 22 8 1 -1.1128 6.6045 2.0284 0.6044 + 23 8 2 0.5564 6.3655 1.0240 0.4799 + 24 8 2 0.5564 5.7494 2.3953 0.7463 + 25 9 1 -1.1128 8.8742 8.0839 4.6354 + 26 9 2 0.5564 9.7892 8.3767 4.4027 + 27 9 2 0.5564 8.3460 8.0813 3.8293 + 28 10 1 -1.1128 5.8097 5.6849 -0.7085 + 29 10 2 0.5564 5.9595 5.0452 -1.5572 + 30 10 2 0.5564 5.2079 6.4617 -0.9303 + 31 11 1 -1.1128 7.6988 2.7233 2.9973 + 32 11 2 0.5564 6.8675 3.0722 3.5888 + 33 11 2 0.5564 7.3299 2.4174 1.9170 + 34 12 1 -1.1128 11.0654 4.5876 3.5549 + 35 12 2 0.5564 10.3773 4.6342 2.8632 + 36 12 2 0.5564 10.4763 5.2035 4.2540 + 37 13 1 -1.1128 2.2356 2.9662 -0.3426 + 38 13 2 0.5564 2.8229 2.0247 -0.6116 + 39 13 2 0.5564 2.9868 3.0933 0.5014 + 40 14 1 -1.1128 0.6871 7.8853 -4.4203 + 41 14 2 0.5564 0.9105 7.8501 -3.3605 + 42 14 2 0.5564 0.0012 8.5407 -4.3902 + 43 15 1 -1.1128 6.1128 -0.4547 -7.5510 + 44 15 2 0.5564 6.0424 -0.9370 -6.6619 + 45 15 2 0.5564 6.9367 -0.2264 -7.6923 + 46 16 1 -1.1128 0.0269 7.9881 0.2323 + 47 16 2 0.5564 0.1935 7.5172 1.0657 + 48 16 2 0.5564 -1.0149 7.8781 0.1064 + 49 17 1 -1.1128 9.0979 4.4375 1.6117 + 50 17 2 0.5564 8.4785 3.5461 2.0370 + 51 17 2 0.5564 8.4603 5.1558 1.1598 + 52 18 1 -1.1128 8.0205 9.5154 -5.6604 + 53 18 2 0.5564 7.0891 9.9739 -5.8289 + 54 18 2 0.5564 8.2865 9.0556 -6.5545 + 55 19 1 -1.1128 7.1473 7.4240 -4.1023 + 56 19 2 0.5564 6.3979 6.9506 -4.4627 + 57 19 2 0.5564 7.1474 8.2756 -4.5903 + 58 20 1 -1.1128 4.2167 1.7495 -7.8409 + 59 20 2 0.5564 5.0855 1.2633 -7.6234 + 60 20 2 0.5564 3.7241 1.1296 -8.5781 + 61 21 1 -1.1128 5.4446 3.8853 -8.9068 + 62 21 2 0.5564 5.1557 4.7500 -8.6367 + 63 21 2 0.5564 4.8182 3.2069 -8.3794 + 64 22 1 -1.1128 6.0386 9.3204 -2.4511 + 65 22 2 0.5564 6.4152 8.5164 -2.8341 + 66 22 2 0.5564 6.7845 9.8389 -1.9444 + 67 23 1 -1.1128 0.1381 2.6793 4.5680 + 68 23 2 0.5564 -0.2809 3.4724 3.9858 + 69 23 2 0.5564 0.2112 2.0940 3.8931 + 70 24 1 -1.1128 11.3965 -3.5321 3.9345 + 71 24 2 0.5564 11.7433 -3.3427 4.8119 + 72 24 2 0.5564 11.9756 -4.4481 3.3455 + 73 25 1 -1.1128 5.1961 -2.3263 -5.0160 + 74 25 2 0.5564 5.5845 -2.6326 -3.9068 + 75 25 2 0.5564 4.9080 -1.3876 -5.0752 + 76 26 1 -1.1128 1.9884 4.6024 -2.2437 + 77 26 2 0.5564 1.1159 4.2702 -2.6316 + 78 26 2 0.5564 2.1782 3.8946 -1.7102 + 79 27 1 -1.1128 3.4286 -3.8254 2.9948 + 80 27 2 0.5564 4.1401 -3.1922 2.7029 + 81 27 2 0.5564 2.6837 -3.3946 2.5442 + 82 28 1 -1.1128 3.4972 4.0312 -4.2534 + 83 28 2 0.5564 3.2969 3.6780 -5.1881 + 84 28 2 0.5564 2.8283 3.9277 -3.5286 + 85 29 1 -1.1128 2.5581 3.6391 -6.8719 + 86 29 2 0.5564 3.0211 2.9680 -7.2340 + 87 29 2 0.5564 1.6628 3.1519 -7.0428 + 88 30 1 -1.1128 5.1229 6.3022 -5.5820 + 89 30 2 0.5564 4.4630 7.1815 -5.7854 + 90 30 2 0.5564 4.3184 5.7177 -5.1800 + 91 31 1 -1.1128 7.5284 -1.4968 -1.3872 + 92 31 2 0.5564 7.6646 -0.7307 -2.1238 + 93 31 2 0.5564 8.3855 -1.4659 -1.0203 + 94 32 1 -1.1128 8.0180 3.6712 -5.5542 + 95 32 2 0.5564 8.2173 4.4555 -6.0837 + 96 32 2 0.5564 8.6346 3.0360 -5.6917 + 97 33 1 -1.1128 16.0706 -0.2590 -5.0003 + 98 33 2 0.5564 16.0528 -0.2534 -5.8842 + 99 33 2 0.5564 15.1342 -0.1953 -4.7488 + 100 34 1 -1.1128 8.7867 0.2163 4.2394 + 101 34 2 0.5564 8.4837 1.0587 3.5818 + 102 34 2 0.5564 9.5314 -0.1659 3.6191 + 103 35 1 -1.1128 3.2350 8.1868 -5.7667 + 104 35 2 0.5564 2.4833 8.4024 -5.2868 + 105 35 2 0.5564 3.9304 8.8555 -5.5968 + 106 36 1 -1.1128 9.4233 5.7040 5.6841 + 107 36 2 0.5564 10.0973 5.7599 6.4736 + 108 36 2 0.5564 8.9638 6.5622 5.3889 + 109 37 1 -1.1128 5.4246 -2.0137 2.6572 + 110 37 2 0.5564 5.8998 -1.6969 3.5099 + 111 37 2 0.5564 5.5652 -1.5574 1.8565 + 112 38 1 -1.1128 -0.1317 1.9336 0.3190 + 113 38 2 0.5564 0.7894 2.2269 -0.1673 + 114 38 2 0.5564 -0.9403 2.4988 0.0245 + 115 39 1 -1.1128 11.8227 3.8276 -2.8190 + 116 39 2 0.5564 11.9045 2.9412 -3.1683 + 117 39 2 0.5564 11.4467 3.8091 -2.0650 + 118 40 1 -1.1128 -1.2368 9.7290 -3.9382 + 119 40 2 0.5564 -1.9257 9.0012 -3.6423 + 120 40 2 0.5564 -1.6406 10.5300 -4.3785 + 121 41 1 -1.1128 10.0754 7.7161 -0.3734 + 122 41 2 0.5564 10.0574 7.7060 -1.5087 + 123 41 2 0.5564 9.2537 7.2615 -0.1815 + 124 42 1 -1.1128 -1.0123 6.0538 -4.7618 + 125 42 2 0.5564 -0.8778 5.3177 -4.1260 + 126 42 2 0.5564 -0.2149 6.6635 -4.6539 + 127 43 1 -1.1128 -1.7545 11.2259 2.6532 + 128 43 2 0.5564 -1.2454 10.3615 3.0272 + 129 43 2 0.5564 -2.3173 10.8721 1.8706 + 130 44 1 -1.1128 12.4899 11.0304 5.7156 + 131 44 2 0.5564 11.9240 11.7099 5.4868 + 132 44 2 0.5564 13.3664 11.1913 5.4911 + 133 45 1 -1.1128 9.6108 7.3259 -3.0932 + 134 45 2 0.5564 8.6400 7.3761 -3.4426 + 135 45 2 0.5564 10.1122 6.7757 -3.8567 + 136 46 1 -1.1128 2.6948 11.0492 4.8330 + 137 46 2 0.5564 2.7024 11.5790 4.0320 + 138 46 2 0.5564 3.0263 9.9915 4.5693 + 139 47 1 -1.1128 2.3981 10.4891 -1.9366 + 140 47 2 0.5564 3.2686 9.9508 -1.8156 + 141 47 2 0.5564 2.0376 10.5907 -0.9310 + 142 48 1 -1.1128 0.9756 10.6255 0.5327 + 143 48 2 0.5564 0.1537 11.3288 0.6051 + 144 48 2 0.5564 0.5959 9.7905 0.6015 + 145 49 1 -1.1128 6.8518 4.4218 -2.9838 + 146 49 2 0.5564 7.5360 4.5692 -4.0439 + 147 49 2 0.5564 6.2902 3.8048 -3.2906 + 148 50 1 -1.1128 5.2164 6.6410 4.0581 + 149 50 2 0.5564 5.2488 6.3928 4.9660 + 150 50 2 0.5564 4.2955 7.0350 3.8882 + 151 51 1 -1.1128 0.3549 6.8498 2.9794 + 152 51 2 0.5564 1.0995 6.6230 3.4814 + 153 51 2 0.5564 -0.3406 6.0721 3.1076 + 154 52 1 -1.1128 8.2297 2.5219 -1.4475 + 155 52 2 0.5564 7.9524 3.3343 -1.9877 + 156 52 2 0.5564 7.7019 2.4769 -0.5961 + 157 53 1 -1.1128 7.8245 6.7224 0.7167 + 158 53 2 0.5564 7.2205 6.3967 0.2032 + 159 53 2 0.5564 7.5432 7.1690 1.5959 + 160 54 1 -1.1128 10.4529 3.5595 -0.4618 + 161 54 2 0.5564 9.6554 3.0113 -1.0854 + 162 54 2 0.5564 9.9860 3.8155 0.4575 + 163 55 1 -1.1128 1.6091 7.3889 -1.9372 + 164 55 2 0.5564 0.8772 7.6442 -1.3852 + 165 55 2 0.5564 1.6148 6.3199 -2.0961 + 166 56 1 -1.1128 4.5847 3.7052 0.8854 + 167 56 2 0.5564 4.9276 4.3327 0.4075 + 168 56 2 0.5564 4.7321 3.6932 1.9221 + 169 57 1 -1.1128 7.7489 0.5850 -3.1580 + 170 57 2 0.5564 8.1265 1.3116 -2.5648 + 171 57 2 0.5564 6.8161 0.9453 -3.2322 + 172 58 1 -1.1128 1.0484 11.8440 -4.0472 + 173 58 2 0.5564 0.7280 11.1323 -4.6797 + 174 58 2 0.5564 1.1330 11.2806 -3.2996 + 175 59 1 -1.1128 4.8434 1.5791 -3.4338 + 176 59 2 0.5564 4.6193 2.5111 -3.8020 + 177 59 2 0.5564 4.5590 0.9918 -4.2585 + 178 60 1 -1.1128 9.5745 -0.8244 7.3550 + 179 60 2 0.5564 9.1589 -1.6054 6.9407 + 180 60 2 0.5564 8.9480 -0.4525 7.9368 + 181 61 1 -1.1128 10.2730 1.8899 6.5393 + 182 61 2 0.5564 10.8570 2.0184 5.8089 + 183 61 2 0.5564 9.9510 0.9418 6.2056 + 184 62 1 -1.1128 4.0753 7.8192 -0.7739 + 185 62 2 0.5564 4.6355 8.0586 -1.5322 + 186 62 2 0.5564 3.2914 7.7124 -1.2255 + 187 63 1 -1.1128 3.7702 0.4562 -1.1286 + 188 63 2 0.5564 3.4391 -0.4159 -1.3938 + 189 63 2 0.5564 4.1240 0.8020 -1.9818 + 190 64 1 -1.1128 6.9868 8.1387 2.7851 + 191 64 2 0.5564 6.2262 7.6175 3.5031 + 192 64 2 0.5564 6.8162 9.0302 2.5521 + +Bonds + + 1 1 1 2 + 2 1 1 3 + 3 1 4 5 + 4 1 4 6 + 5 1 7 8 + 6 1 7 9 + 7 1 10 11 + 8 1 10 12 + 9 1 13 14 + 10 1 13 15 + 11 1 16 17 + 12 1 16 18 + 13 1 19 20 + 14 1 19 21 + 15 1 22 23 + 16 1 22 24 + 17 1 25 26 + 18 1 25 27 + 19 1 28 29 + 20 1 28 30 + 21 1 31 32 + 22 1 31 33 + 23 1 34 35 + 24 1 34 36 + 25 1 37 38 + 26 1 37 39 + 27 1 40 41 + 28 1 40 42 + 29 1 43 44 + 30 1 43 45 + 31 1 46 47 + 32 1 46 48 + 33 1 49 50 + 34 1 49 51 + 35 1 52 53 + 36 1 52 54 + 37 1 55 56 + 38 1 55 57 + 39 1 58 59 + 40 1 58 60 + 41 1 61 62 + 42 1 61 63 + 43 1 64 65 + 44 1 64 66 + 45 1 67 68 + 46 1 67 69 + 47 1 70 71 + 48 1 70 72 + 49 1 73 74 + 50 1 73 75 + 51 1 76 77 + 52 1 76 78 + 53 1 79 80 + 54 1 79 81 + 55 1 82 83 + 56 1 82 84 + 57 1 85 86 + 58 1 85 87 + 59 1 88 89 + 60 1 88 90 + 61 1 91 92 + 62 1 91 93 + 63 1 94 95 + 64 1 94 96 + 65 1 97 98 + 66 1 97 99 + 67 1 100 101 + 68 1 100 102 + 69 1 103 104 + 70 1 103 105 + 71 1 106 107 + 72 1 106 108 + 73 1 109 110 + 74 1 109 111 + 75 1 112 113 + 76 1 112 114 + 77 1 115 116 + 78 1 115 117 + 79 1 118 119 + 80 1 118 120 + 81 1 121 122 + 82 1 121 123 + 83 1 124 125 + 84 1 124 126 + 85 1 127 128 + 86 1 127 129 + 87 1 130 131 + 88 1 130 132 + 89 1 133 134 + 90 1 133 135 + 91 1 136 137 + 92 1 136 138 + 93 1 139 140 + 94 1 139 141 + 95 1 142 143 + 96 1 142 144 + 97 1 145 146 + 98 1 145 147 + 99 1 148 149 + 100 1 148 150 + 101 1 151 152 + 102 1 151 153 + 103 1 154 155 + 104 1 154 156 + 105 1 157 158 + 106 1 157 159 + 107 1 160 161 + 108 1 160 162 + 109 1 163 164 + 110 1 163 165 + 111 1 166 167 + 112 1 166 168 + 113 1 169 170 + 114 1 169 171 + 115 1 172 173 + 116 1 172 174 + 117 1 175 176 + 118 1 175 177 + 119 1 178 179 + 120 1 178 180 + 121 1 181 182 + 122 1 181 183 + 123 1 184 185 + 124 1 184 186 + 125 1 187 188 + 126 1 187 189 + 127 1 190 191 + 128 1 190 192 + +Angles + + 1 1 2 1 3 + 2 1 5 4 6 + 3 1 8 7 9 + 4 1 11 10 12 + 5 1 14 13 15 + 6 1 17 16 18 + 7 1 20 19 21 + 8 1 23 22 24 + 9 1 26 25 27 + 10 1 29 28 30 + 11 1 32 31 33 + 12 1 35 34 36 + 13 1 38 37 39 + 14 1 41 40 42 + 15 1 44 43 45 + 16 1 47 46 48 + 17 1 50 49 51 + 18 1 53 52 54 + 19 1 56 55 57 + 20 1 59 58 60 + 21 1 62 61 63 + 22 1 65 64 66 + 23 1 68 67 69 + 24 1 71 70 72 + 25 1 74 73 75 + 26 1 77 76 78 + 27 1 80 79 81 + 28 1 83 82 84 + 29 1 86 85 87 + 30 1 89 88 90 + 31 1 92 91 93 + 32 1 95 94 96 + 33 1 98 97 99 + 34 1 101 100 102 + 35 1 104 103 105 + 36 1 107 106 108 + 37 1 110 109 111 + 38 1 113 112 114 + 39 1 116 115 117 + 40 1 119 118 120 + 41 1 122 121 123 + 42 1 125 124 126 + 43 1 128 127 129 + 44 1 131 130 132 + 45 1 134 133 135 + 46 1 137 136 138 + 47 1 140 139 141 + 48 1 143 142 144 + 49 1 146 145 147 + 50 1 149 148 150 + 51 1 152 151 153 + 52 1 155 154 156 + 53 1 158 157 159 + 54 1 161 160 162 + 55 1 164 163 165 + 56 1 167 166 168 + 57 1 170 169 171 + 58 1 173 172 174 + 59 1 176 175 177 + 60 1 179 178 180 + 61 1 182 181 183 + 62 1 185 184 186 + 63 1 188 187 189 + 64 1 191 190 192 + diff --git a/examples/USER/gle/in.tip4pf b/examples/USER/gle/in.tip4pf new file mode 100644 index 0000000000..606042f35a --- /dev/null +++ b/examples/USER/gle/in.tip4pf @@ -0,0 +1,29 @@ +units real +atom_style gle # 6 + +pair_style lj/cut/tip4p/long 1 2 1 1 0.14714951 8 +bond_style class2 +angle_style harmonic +kspace_style pppm 0.0001 + +read_data data.tip4pf +pair_coeff * 2 0 0 +pair_coeff 1 1 0.1852 3.1589022 + +#q-tip4pf bond parameters +bond_coeff 1 0.9419 607.19354 -1388.6516 1852.577 +angle_coeff 1 43.93 107.4 +thermo 100 +thermo_style custom step temp pe ke etotal + +timestep 0.5 + +dump snapxyz all atom 100 h2o.lammpstraj + +# some problem +fix 1 all gle 6 smart.A 300 31415 +fix_modify 1 energy yes +#fix 1 all nve + +run 10000 + diff --git a/src/Makefile b/src/Makefile index f8e70a94dc..445045914f 100644 --- a/src/Makefile +++ b/src/Makefile @@ -19,7 +19,7 @@ PACKAGE = asphere body class2 colloid dipole fld gpu granular kim \ PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \ user-cuda user-eff user-fep user-lb user-misc user-molfile \ - user-omp user-phonon user-qmmm user-reaxc user-sph + user-omp user-phonon user-qmmm user-reaxc user-sph user-gle PACKLIB = gpu kim meam poems reax voronoi \ user-atc user-awpmd user-colvars user-qmmm user-cuda user-molfile diff --git a/src/USER-GLE/atom_vec_gle.cpp b/src/USER-GLE/atom_vec_gle.cpp new file mode 100644 index 0000000000..c83e2c6ea6 --- /dev/null +++ b/src/USER-GLE/atom_vec_gle.cpp @@ -0,0 +1,1161 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "atom_vec_gle.h" +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "modify.h" +#include "fix.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* GLE momenta should be carried around with the atoms as they move across + * processors. Don't know what is the correct place to do that.... + */ +#define COMM_GLES // communicates GLE additional momenta across processors +#define DELTA 10000 + +/* ---------------------------------------------------------------------- */ + +AtomVecGLE::AtomVecGLE(LAMMPS * lmp): AtomVec(lmp) //, int narg, char **arg) : AtomVec(lmp, narg, arg) +{ + + //! AtomVec classes used to have the infrastructure for reading extra arguments. Now I am hardcoding this, but it should really be a parameter. + //if (narg>0) + { + ns = 6; //atoi(arg[0]); + } + atom->ns=ns; + + molecular = 1; + bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 1; + mass_type = 1; + + comm_x_only = comm_f_only = 1; + size_forward = 3; + size_reverse = 3; + size_border = 8; +#ifdef COMM_GLES + size_velocity = 3 + 3 * (atom->ns+1); // also transfer additional momenta with velocities! +#else + size_velocity = 3; +#endif + size_data_atom = 7; + size_data_vel = 4; // may need to manipulate this in order to save additional momenta to data file + xcol_data = 5; + + atom->molecule_flag = atom->q_flag = 1; +} + +/* ---------------------------------------------------------------------- + grow atom arrays + n = 0 grows arrays by DELTA + n > 0 allocates arrays to size n +------------------------------------------------------------------------- */ + +void AtomVecGLE::grow(int n) +{ + if (n == 0) nmax += DELTA; + else nmax = n; + atom->nmax = nmax; + if (nmax < 0 || nmax > MAXSMALLINT) + error->one(FLERR,"Per-processor system is too big"); + + tag = memory->grow(atom->tag,nmax,"atom:tag"); + type = memory->grow(atom->type,nmax,"atom:type"); + mask = memory->grow(atom->mask,nmax,"atom:mask"); + image = memory->grow(atom->image,nmax,"atom:image"); + x = memory->grow(atom->x,nmax,3,"atom:x"); + v = memory->grow(atom->v,nmax,3,"atom:v"); + f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f"); + s = memory->grow(atom->s,nmax*comm->nthreads,3*(ns+1),"atom:s"); + + q = memory->grow(atom->q,nmax,"atom:q"); + molecule = memory->grow(atom->molecule,nmax,"atom:molecule"); + + nspecial = memory->grow(atom->nspecial,nmax,3,"atom:nspecial"); + special = memory->grow(atom->special,nmax,atom->maxspecial,"atom:special"); + + num_bond = memory->grow(atom->num_bond,nmax,"atom:num_bond"); + bond_type = memory->grow(atom->bond_type,nmax,atom->bond_per_atom, + "atom:bond_type"); + bond_atom = memory->grow(atom->bond_atom,nmax,atom->bond_per_atom, + "atom:bond_atom"); + + num_angle = memory->grow(atom->num_angle,nmax,"atom:num_angle"); + angle_type = memory->grow(atom->angle_type,nmax,atom->angle_per_atom, + "atom:angle_type"); + angle_atom1 = memory->grow(atom->angle_atom1,nmax,atom->angle_per_atom, + "atom:angle_atom1"); + angle_atom2 = memory->grow(atom->angle_atom2,nmax,atom->angle_per_atom, + "atom:angle_atom2"); + angle_atom3 = memory->grow(atom->angle_atom3,nmax,atom->angle_per_atom, + "atom:angle_atom3"); + + num_dihedral = memory->grow(atom->num_dihedral,nmax,"atom:num_dihedral"); + dihedral_type = memory->grow(atom->dihedral_type,nmax, + atom->dihedral_per_atom,"atom:dihedral_type"); + dihedral_atom1 = + memory->grow(atom->dihedral_atom1,nmax,atom->dihedral_per_atom, + "atom:dihedral_atom1"); + dihedral_atom2 = + memory->grow(atom->dihedral_atom2,nmax,atom->dihedral_per_atom, + "atom:dihedral_atom2"); + dihedral_atom3 = + memory->grow(atom->dihedral_atom3,nmax,atom->dihedral_per_atom, + "atom:dihedral_atom3"); + dihedral_atom4 = + memory->grow(atom->dihedral_atom4,nmax,atom->dihedral_per_atom, + "atom:dihedral_atom4"); + + num_improper = memory->grow(atom->num_improper,nmax,"atom:num_improper"); + improper_type = + memory->grow(atom->improper_type,nmax,atom->improper_per_atom, + "atom:improper_type"); + improper_atom1 = + memory->grow(atom->improper_atom1,nmax,atom->improper_per_atom, + "atom:improper_atom1"); + improper_atom2 = + memory->grow(atom->improper_atom2,nmax,atom->improper_per_atom, + "atom:improper_atom2"); + improper_atom3 = + memory->grow(atom->improper_atom3,nmax,atom->improper_per_atom, + "atom:improper_atom3"); + improper_atom4 = + memory->grow(atom->improper_atom4,nmax,atom->improper_per_atom, + "atom:improper_atom4"); + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); +} + +/* ---------------------------------------------------------------------- + reset local array ptrs +------------------------------------------------------------------------- */ + +void AtomVecGLE::grow_reset() +{ + tag = atom->tag; type = atom->type; + mask = atom->mask; image = atom->image; + x = atom->x; v = atom->v; f = atom->f; s = atom->s; + q = atom->q; molecule = atom->molecule; + nspecial = atom->nspecial; special = atom->special; + num_bond = atom->num_bond; bond_type = atom->bond_type; + bond_atom = atom->bond_atom; + num_angle = atom->num_angle; angle_type = atom->angle_type; + angle_atom1 = atom->angle_atom1; angle_atom2 = atom->angle_atom2; + angle_atom3 = atom->angle_atom3; + num_dihedral = atom->num_dihedral; dihedral_type = atom->dihedral_type; + dihedral_atom1 = atom->dihedral_atom1; dihedral_atom2 = atom->dihedral_atom2; + dihedral_atom3 = atom->dihedral_atom3; dihedral_atom4 = atom->dihedral_atom4; + num_improper = atom->num_improper; improper_type = atom->improper_type; + improper_atom1 = atom->improper_atom1; improper_atom2 = atom->improper_atom2; + improper_atom3 = atom->improper_atom3; improper_atom4 = atom->improper_atom4; +} + +/* ---------------------------------------------------------------------- + copy atom I info to atom J +------------------------------------------------------------------------- */ + +void AtomVecGLE::copy(int i, int j, int delflag) +{ + int k; + + tag[j] = tag[i]; + type[j] = type[i]; + mask[j] = mask[i]; + image[j] = image[i]; + x[j][0] = x[i][0]; + x[j][1] = x[i][1]; + x[j][2] = x[i][2]; + v[j][0] = v[i][0]; + v[j][1] = v[i][1]; + v[j][2] = v[i][2]; + + for (k = 0; k < 3*(ns+1); k++) s[j][k] = s[i][k]; + + q[j] = q[i]; + molecule[j] = molecule[i]; + + num_bond[j] = num_bond[i]; + for (k = 0; k < num_bond[j]; k++) { + bond_type[j][k] = bond_type[i][k]; + bond_atom[j][k] = bond_atom[i][k]; + } + + num_angle[j] = num_angle[i]; + for (k = 0; k < num_angle[j]; k++) { + angle_type[j][k] = angle_type[i][k]; + angle_atom1[j][k] = angle_atom1[i][k]; + angle_atom2[j][k] = angle_atom2[i][k]; + angle_atom3[j][k] = angle_atom3[i][k]; + } + + num_dihedral[j] = num_dihedral[i]; + for (k = 0; k < num_dihedral[j]; k++) { + dihedral_type[j][k] = dihedral_type[i][k]; + dihedral_atom1[j][k] = dihedral_atom1[i][k]; + dihedral_atom2[j][k] = dihedral_atom2[i][k]; + dihedral_atom3[j][k] = dihedral_atom3[i][k]; + dihedral_atom4[j][k] = dihedral_atom4[i][k]; + } + + num_improper[j] = num_improper[i]; + for (k = 0; k < num_improper[j]; k++) { + improper_type[j][k] = improper_type[i][k]; + improper_atom1[j][k] = improper_atom1[i][k]; + improper_atom2[j][k] = improper_atom2[i][k]; + improper_atom3[j][k] = improper_atom3[i][k]; + improper_atom4[j][k] = improper_atom4[i][k]; + } + + nspecial[j][0] = nspecial[i][0]; + nspecial[j][1] = nspecial[i][1]; + nspecial[j][2] = nspecial[i][2]; + for (k = 0; k < nspecial[j][2]; k++) special[j][k] = special[i][k]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j,delflag); +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecGLE::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecGLE::pack_comm_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,k,m; + double dx,dy,dz,dvx,dvy,dvz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; +#ifdef COMM_GLES + for (k=0; k<3*(ns+1); k++) buf[m++] = s[j][k]; +#endif + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; + dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; + dz = pbc[2]*domain->zprd; + } + if (!deform_vremap) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; +#ifdef COMM_GLES + for (k=0; k<3*(ns+1); k++) buf[m++] = s[j][k]; +#endif + } + } else { + dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4]; + dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3]; + dvz = pbc[2]*h_rate[2]; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + if (mask[i] & deform_groupbit) { + buf[m++] = v[j][0] + dvx; + buf[m++] = v[j][1] + dvy; + buf[m++] = v[j][2] + dvz; + } else { + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } +#ifdef COMM_GLES + for (k=0; k<3*(ns+1); k++) buf[m++] = s[j][k]; +#endif + } + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecGLE::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecGLE::unpack_comm_vel(int n, int first, double *buf) +{ + int i,k,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; +#ifdef COMM_GLES + for (k=0; k<3*(ns+1); k++) s[i][k] = buf[m++]; +#endif + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecGLE::pack_reverse(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = f[i][0]; + buf[m++] = f[i][1]; + buf[m++] = f[i][2]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecGLE::unpack_reverse(int n, int *list, double *buf) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + f[j][0] += buf[m++]; + f[j][1] += buf[m++]; + f[j][2] += buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecGLE::pack_border(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = ubuf(tag[j]).d; + buf[m++] = ubuf(type[j]).d; + buf[m++] = ubuf(mask[j]).d; + buf[m++] = q[j]; + buf[m++] = ubuf(molecule[j]).d; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = ubuf(tag[j]).d; + buf[m++] = ubuf(type[j]).d; + buf[m++] = ubuf(mask[j]).d; + buf[m++] = q[j]; + buf[m++] = ubuf(molecule[j]).d; + } + } + + if (atom->nextra_border) + for (int iextra = 0; iextra < atom->nextra_border; iextra++) + m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]); + + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecGLE::pack_border_vel(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,k,m; + double dx,dy,dz,dvx,dvy,dvz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = ubuf(tag[j]).d; + buf[m++] = ubuf(type[j]).d; + buf[m++] = ubuf(mask[j]).d; + buf[m++] = q[j]; + buf[m++] = ubuf(molecule[j]).d; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; +#ifdef COMM_GLES + for (k=0; k<3*(ns+1); k++) buf[m++] = s[j][k]; +#endif + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + if (!deform_vremap) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = ubuf(tag[j]).d; + buf[m++] = ubuf(type[j]).d; + buf[m++] = ubuf(mask[j]).d; + buf[m++] = q[j]; + buf[m++] = ubuf(molecule[j]).d; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; +#ifdef COMM_GLES + for (k=0; k<3*(ns+1); k++) buf[m++] = s[j][k]; +#endif + } + } else { + dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4]; + dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3]; + dvz = pbc[2]*h_rate[2]; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = ubuf(tag[j]).d; + buf[m++] = ubuf(type[j]).d; + buf[m++] = ubuf(mask[j]).d; + buf[m++] = q[j]; + buf[m++] = ubuf(molecule[j]).d; + if (mask[i] & deform_groupbit) { + buf[m++] = v[j][0] + dvx; + buf[m++] = v[j][1] + dvy; + buf[m++] = v[j][2] + dvz; + } else { + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + } +#ifdef COMM_GLES + for (k=0; k<3*(ns+1); k++) buf[m++] = s[j][k]; +#endif + } + } + } + + if (atom->nextra_border) + for (int iextra = 0; iextra < atom->nextra_border; iextra++) + m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]); + + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecGLE::pack_border_hybrid(int n, int *list, double *buf) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = q[j]; + buf[m++] = ubuf(molecule[j]).d; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecGLE::unpack_border(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = (tagint) ubuf(buf[m++]).i; + type[i] = (int) ubuf(buf[m++]).i; + mask[i] = (int) ubuf(buf[m++]).i; + q[i] = buf[m++]; + molecule[i] = (tagint) ubuf(buf[m++]).i; + } + + if (atom->nextra_border) + for (int iextra = 0; iextra < atom->nextra_border; iextra++) + m += modify->fix[atom->extra_border[iextra]]-> + unpack_border(n,first,&buf[m]); +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecGLE::unpack_border_vel(int n, int first, double *buf) +{ + int i,k,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = (tagint) ubuf(buf[m++]).i; + type[i] = (int) ubuf(buf[m++]).i; + mask[i] = (int) ubuf(buf[m++]).i; + q[i] = buf[m++]; + molecule[i] = (tagint) ubuf(buf[m++]).i; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; +#ifdef COMM_GLES + for (k=0; k<3*(ns+1); k++) s[i][k] = buf[m++]; +#endif + } + + if (atom->nextra_border) + for (int iextra = 0; iextra < atom->nextra_border; iextra++) + m += modify->fix[atom->extra_border[iextra]]-> + unpack_border(n,first,&buf[m]); +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecGLE::unpack_border_hybrid(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + q[i] = buf[m++]; + molecule[i] = (tagint) ubuf(buf[m++]).i; + } + return m; +} + +/* ---------------------------------------------------------------------- + pack data for atom I for sending to another proc + xyz must be 1st 3 values, so comm::exchange() can test on them +------------------------------------------------------------------------- */ + +int AtomVecGLE::pack_exchange(int i, double *buf) +{ + int k; + + int m = 1; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + buf[m++] = v[i][0]; + buf[m++] = v[i][1]; + buf[m++] = v[i][2]; + buf[m++] = ubuf(tag[i]).d; + buf[m++] = ubuf(type[i]).d; + buf[m++] = ubuf(mask[i]).d; + buf[m++] = ubuf(image[i]).d; + + buf[m++] = q[i]; + buf[m++] = ubuf(molecule[i]).d; + + buf[m++] = ubuf(num_bond[i]).d; + for (k = 0; k < num_bond[i]; k++) { + buf[m++] = ubuf(bond_type[i][k]).d; + buf[m++] = ubuf(bond_atom[i][k]).d; + } + + buf[m++] = ubuf(num_angle[i]).d; + for (k = 0; k < num_angle[i]; k++) { + buf[m++] = ubuf(angle_type[i][k]).d; + buf[m++] = ubuf(angle_atom1[i][k]).d; + buf[m++] = ubuf(angle_atom2[i][k]).d; + buf[m++] = ubuf(angle_atom3[i][k]).d; + } + + buf[m++] = ubuf(num_dihedral[i]).d; + for (k = 0; k < num_dihedral[i]; k++) { + buf[m++] = ubuf(dihedral_type[i][k]).d; + buf[m++] = ubuf(dihedral_atom1[i][k]).d; + buf[m++] = ubuf(dihedral_atom2[i][k]).d; + buf[m++] = ubuf(dihedral_atom3[i][k]).d; + buf[m++] = ubuf(dihedral_atom4[i][k]).d; + } + + buf[m++] = ubuf(num_improper[i]).d; + for (k = 0; k < num_improper[i]; k++) { + buf[m++] = ubuf(improper_type[i][k]).d; + buf[m++] = ubuf(improper_atom1[i][k]).d; + buf[m++] = ubuf(improper_atom2[i][k]).d; + buf[m++] = ubuf(improper_atom3[i][k]).d; + buf[m++] = ubuf(improper_atom4[i][k]).d; + } + + buf[m++] = ubuf(nspecial[i][0]).d; + buf[m++] = ubuf(nspecial[i][1]).d; + buf[m++] = ubuf(nspecial[i][2]).d; + for (k = 0; k < nspecial[i][2]; k++) buf[m++] = ubuf(special[i][k]).d; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]); + + buf[0] = m; + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecGLE::unpack_exchange(double *buf) +{ + int k; + + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + int m = 1; + x[nlocal][0] = buf[m++]; + x[nlocal][1] = buf[m++]; + x[nlocal][2] = buf[m++]; + v[nlocal][0] = buf[m++]; + v[nlocal][1] = buf[m++]; + v[nlocal][2] = buf[m++]; +#ifdef COMM_GLES + for (k=0; k<3*(ns+1); k++) s[nlocal][k] = buf[m++]; +#endif + + tag[nlocal] = (tagint) ubuf(buf[m++]).i; + type[nlocal] = (int) ubuf(buf[m++]).i; + mask[nlocal] = (int) ubuf(buf[m++]).i; + image[nlocal] = (imageint) ubuf(buf[m++]).i; + + q[nlocal] = buf[m++]; + molecule[nlocal] = (tagint) ubuf(buf[m++]).i; + + num_bond[nlocal] = (int) ubuf(buf[m++]).i; + for (k = 0; k < num_bond[nlocal]; k++) { + bond_type[nlocal][k] = (int) ubuf(buf[m++]).i; + bond_atom[nlocal][k] = (tagint) ubuf(buf[m++]).i; + } + + num_angle[nlocal] = (int) ubuf(buf[m++]).i; + for (k = 0; k < num_angle[nlocal]; k++) { + angle_type[nlocal][k] = (int) ubuf(buf[m++]).i; + angle_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i; + angle_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i; + angle_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i; + } + + num_dihedral[nlocal] = (int) ubuf(buf[m++]).i; + for (k = 0; k < num_dihedral[nlocal]; k++) { + dihedral_type[nlocal][k] = (int) ubuf(buf[m++]).i; + dihedral_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i; + dihedral_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i; + dihedral_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i; + dihedral_atom4[nlocal][k] = (tagint) ubuf(buf[m++]).i; + } + + num_improper[nlocal] = (int) ubuf(buf[m++]).i; + for (k = 0; k < num_improper[nlocal]; k++) { + improper_type[nlocal][k] = (int) ubuf(buf[m++]).i; + improper_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i; + improper_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i; + improper_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i; + improper_atom4[nlocal][k] = (tagint) ubuf(buf[m++]).i; + } + + nspecial[nlocal][0] = (int) ubuf(buf[m++]).i; + nspecial[nlocal][1] = (int) ubuf(buf[m++]).i; + nspecial[nlocal][2] = (int) ubuf(buf[m++]).i; + for (k = 0; k < nspecial[nlocal][2]; k++) + special[nlocal][k] = (tagint) ubuf(buf[m++]).i; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + m += modify->fix[atom->extra_grow[iextra]]-> + unpack_exchange(nlocal,&buf[m]); + + atom->nlocal++; + return m; +} + +/* ---------------------------------------------------------------------- + size of restart data for all atoms owned by this proc + include extra data stored by fixes +------------------------------------------------------------------------- */ + +int AtomVecGLE::size_restart() +{ + int i; + + int nlocal = atom->nlocal; + int n = 0; + for (i = 0; i < nlocal; i++) + n += 17 + 2*num_bond[i] + 4*num_angle[i] + + 5*num_dihedral[i] + 5*num_improper[i]; + + if (atom->nextra_restart) + for (int iextra = 0; iextra < atom->nextra_restart; iextra++) + for (i = 0; i < nlocal; i++) + n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); + + return n; +} + +/* ---------------------------------------------------------------------- + pack atom I's data for restart file including extra quantities + xyz must be 1st 3 values, so that read_restart can test on them + molecular types may be negative, but write as positive +------------------------------------------------------------------------- */ + +int AtomVecGLE::pack_restart(int i, double *buf) +{ + int k; + + int m = 1; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + buf[m++] = ubuf(tag[i]).d; + buf[m++] = ubuf(type[i]).d; + buf[m++] = ubuf(mask[i]).d; + buf[m++] = ubuf(image[i]).d; + buf[m++] = v[i][0]; + buf[m++] = v[i][1]; + buf[m++] = v[i][2]; + + buf[m++] = q[i]; + buf[m++] = ubuf(molecule[i]).d; + + buf[m++] = ubuf(num_bond[i]).d; + for (k = 0; k < num_bond[i]; k++) { + buf[m++] = ubuf(MAX(bond_type[i][k],-bond_type[i][k])).d; + buf[m++] = ubuf(bond_atom[i][k]).d; + } + + buf[m++] = ubuf(num_angle[i]).d; + for (k = 0; k < num_angle[i]; k++) { + buf[m++] = ubuf(MAX(angle_type[i][k],-angle_type[i][k])).d; + buf[m++] = ubuf(angle_atom1[i][k]).d; + buf[m++] = ubuf(angle_atom2[i][k]).d; + buf[m++] = ubuf(angle_atom3[i][k]).d; + } + + buf[m++] = ubuf(num_dihedral[i]).d; + for (k = 0; k < num_dihedral[i]; k++) { + buf[m++] = ubuf(MAX(dihedral_type[i][k],-dihedral_type[i][k])).d; + buf[m++] = ubuf(dihedral_atom1[i][k]).d; + buf[m++] = ubuf(dihedral_atom2[i][k]).d; + buf[m++] = ubuf(dihedral_atom3[i][k]).d; + buf[m++] = ubuf(dihedral_atom4[i][k]).d; + } + + buf[m++] = ubuf(num_improper[i]).d; + for (k = 0; k < num_improper[i]; k++) { + buf[m++] = ubuf(MAX(improper_type[i][k],-improper_type[i][k])).d; + buf[m++] = ubuf(improper_atom1[i][k]).d; + buf[m++] = ubuf(improper_atom2[i][k]).d; + buf[m++] = ubuf(improper_atom3[i][k]).d; + buf[m++] = ubuf(improper_atom4[i][k]).d; + } + + if (atom->nextra_restart) + for (int iextra = 0; iextra < atom->nextra_restart; iextra++) + m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]); + + buf[0] = m; + return m; +} + +/* ---------------------------------------------------------------------- + unpack data for one atom from restart file including extra quantities +------------------------------------------------------------------------- */ + +int AtomVecGLE::unpack_restart(double *buf) +{ + int k; + + int nlocal = atom->nlocal; + if (nlocal == nmax) { + grow(0); + if (atom->nextra_store) + memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra"); + } + + int m = 1; + x[nlocal][0] = buf[m++]; + x[nlocal][1] = buf[m++]; + x[nlocal][2] = buf[m++]; + tag[nlocal] = (tagint) ubuf(buf[m++]).i; + type[nlocal] = (int) ubuf(buf[m++]).i; + mask[nlocal] = (int) ubuf(buf[m++]).i; + image[nlocal] = (imageint) ubuf(buf[m++]).i; + v[nlocal][0] = buf[m++]; + v[nlocal][1] = buf[m++]; + v[nlocal][2] = buf[m++]; + + q[nlocal] = buf[m++]; + molecule[nlocal] = (tagint) ubuf(buf[m++]).i; + + num_bond[nlocal] = (int) ubuf(buf[m++]).i; + for (k = 0; k < num_bond[nlocal]; k++) { + bond_type[nlocal][k] = (int) ubuf(buf[m++]).i; + bond_atom[nlocal][k] = (tagint) ubuf(buf[m++]).i; + } + + num_angle[nlocal] = (int) ubuf(buf[m++]).i; + for (k = 0; k < num_angle[nlocal]; k++) { + angle_type[nlocal][k] = (int) ubuf(buf[m++]).i; + angle_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i; + angle_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i; + angle_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i; + } + + num_dihedral[nlocal] = (int) ubuf(buf[m++]).i; + for (k = 0; k < num_dihedral[nlocal]; k++) { + dihedral_type[nlocal][k] = (int) ubuf(buf[m++]).i; + dihedral_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i; + dihedral_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i; + dihedral_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i; + dihedral_atom4[nlocal][k] = (tagint) ubuf(buf[m++]).i; + } + + num_improper[nlocal] = (int) ubuf(buf[m++]).i; + for (k = 0; k < num_improper[nlocal]; k++) { + improper_type[nlocal][k] = (int) ubuf(buf[m++]).i; + improper_atom1[nlocal][k] = (tagint) ubuf(buf[m++]).i; + improper_atom2[nlocal][k] = (tagint) ubuf(buf[m++]).i; + improper_atom3[nlocal][k] = (tagint) ubuf(buf[m++]).i; + improper_atom4[nlocal][k] = (tagint) ubuf(buf[m++]).i; + } + + nspecial[nlocal][0] = nspecial[nlocal][1] = nspecial[nlocal][2] = 0; + + double **extra = atom->extra; + if (atom->nextra_store) { + int size = static_cast (buf[0]) - m; + for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++]; + } + + atom->nlocal++; + return m; +} + +/* ---------------------------------------------------------------------- + create one atom of itype at coord + set other values to defaults +------------------------------------------------------------------------- */ + +void AtomVecGLE::create_atom(int itype, double *coord) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + tag[nlocal] = 0; + type[nlocal] = itype; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + mask[nlocal] = 1; + image[nlocal] = ((imageint) IMGMAX << IMG2BITS) | + ((imageint) IMGMAX << IMGBITS) | IMGMAX; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + for (int k=0; k<3*(ns+1); k++) s[nlocal][k] = 0.0; + + q[nlocal] = 0.0; + molecule[nlocal] = 0; + num_bond[nlocal] = 0; + num_angle[nlocal] = 0; + num_dihedral[nlocal] = 0; + num_improper[nlocal] = 0; + nspecial[nlocal][0] = nspecial[nlocal][1] = nspecial[nlocal][2] = 0; + + atom->nlocal++; +} + +/* ---------------------------------------------------------------------- + unpack one line from Atoms section of data file + initialize other atom quantities +------------------------------------------------------------------------- */ + +void AtomVecGLE::data_atom(double *coord, tagint imagetmp, char **values) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + tag[nlocal] = ATOTAGINT(values[0]); + molecule[nlocal] = ATOTAGINT(values[1]); + type[nlocal] = atoi(values[2]); + if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) + error->one(FLERR,"Invalid atom type in Atoms section of data file"); + + q[nlocal] = atof(values[3]); + + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + + image[nlocal] = imagetmp; + + mask[nlocal] = 1; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + for (int k=0; k<3*(ns+1); k++) s[nlocal][k] = 0.0; + num_bond[nlocal] = 0; + num_angle[nlocal] = 0; + num_dihedral[nlocal] = 0; + num_improper[nlocal] = 0; + + atom->nlocal++; +} + +/* ---------------------------------------------------------------------- + unpack hybrid quantities from one line in Atoms section of data file + initialize other atom quantities for this sub-style +------------------------------------------------------------------------- */ + +int AtomVecGLE::data_atom_hybrid(int nlocal, char **values) +{ + molecule[nlocal] = ATOTAGINT(values[0]); + q[nlocal] = atof(values[1]); + + num_bond[nlocal] = 0; + num_angle[nlocal] = 0; + num_dihedral[nlocal] = 0; + num_improper[nlocal] = 0; + + return 2; +} + +/* ---------------------------------------------------------------------- + pack atom info for data file including 3 image flags +------------------------------------------------------------------------- */ + +void AtomVecGLE::pack_data(double **buf) +{ + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) { + buf[i][0] = ubuf(tag[i]).d; + buf[i][1] = ubuf(molecule[i]).d; + buf[i][2] = ubuf(type[i]).d; + buf[i][3] = q[i]; + buf[i][4] = x[i][0]; + buf[i][5] = x[i][1]; + buf[i][6] = x[i][2]; + buf[i][7] = ubuf((image[i] & IMGMASK) - IMGMAX).d; + buf[i][8] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d; + buf[i][9] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d; + } +} + +/* ---------------------------------------------------------------------- + pack hybrid atom info for data file +------------------------------------------------------------------------- */ + +int AtomVecGLE::pack_data_hybrid(int i, double *buf) +{ + buf[0] = ubuf(molecule[i]).d; + buf[1] = q[i]; + return 2; +} + +/* ---------------------------------------------------------------------- + write atom info to data file including 3 image flags +------------------------------------------------------------------------- */ + +void AtomVecGLE::write_data(FILE *fp, int n, double **buf) +{ + for (int i = 0; i < n; i++) + fprintf(fp,TAGINT_FORMAT " " TAGINT_FORMAT + " %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n", + (tagint) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i, + (int) ubuf(buf[i][2]).i, + buf[i][3],buf[i][4],buf[i][5],buf[i][6], + (int) ubuf(buf[i][7]).i,(int) ubuf(buf[i][8]).i, + (int) ubuf(buf[i][9]).i); +} + +/* ---------------------------------------------------------------------- + write hybrid atom info to data file +------------------------------------------------------------------------- */ + +int AtomVecGLE::write_data_hybrid(FILE *fp, double *buf) +{ + fprintf(fp," " TAGINT_FORMAT " %-1.16e",(tagint) ubuf(buf[0]).i,buf[1]); + return 2; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory +------------------------------------------------------------------------- */ + +bigint AtomVecGLE::memory_usage() +{ + bigint bytes = 0; + + if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax); + if (atom->memcheck("type")) bytes += memory->usage(type,nmax); + if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax); + if (atom->memcheck("image")) bytes += memory->usage(image,nmax); + if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3); + if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3); + if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3); + if (atom->memcheck("s")) bytes += memory->usage(f,nmax,3*(ns+1)); + + + if (atom->memcheck("q")) bytes += memory->usage(q,nmax); + if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax); + if (atom->memcheck("nspecial")) bytes += memory->usage(nspecial,nmax,3); + if (atom->memcheck("special")) + bytes += memory->usage(special,nmax,atom->maxspecial); + + if (atom->memcheck("num_bond")) bytes += memory->usage(num_bond,nmax); + if (atom->memcheck("bond_type")) + bytes += memory->usage(bond_type,nmax,atom->bond_per_atom); + if (atom->memcheck("bond_atom")) + bytes += memory->usage(bond_atom,nmax,atom->bond_per_atom); + + if (atom->memcheck("num_angle")) bytes += memory->usage(num_angle,nmax); + if (atom->memcheck("angle_type")) + bytes += memory->usage(angle_type,nmax,atom->angle_per_atom); + if (atom->memcheck("angle_atom1")) + bytes += memory->usage(angle_atom1,nmax,atom->angle_per_atom); + if (atom->memcheck("angle_atom2")) + bytes += memory->usage(angle_atom2,nmax,atom->angle_per_atom); + if (atom->memcheck("angle_atom3")) + bytes += memory->usage(angle_atom3,nmax,atom->angle_per_atom); + + if (atom->memcheck("num_dihedral")) bytes += memory->usage(num_dihedral,nmax); + if (atom->memcheck("dihedral_type")) + bytes += memory->usage(dihedral_type,nmax,atom->dihedral_per_atom); + if (atom->memcheck("dihedral_atom1")) + bytes += memory->usage(dihedral_atom1,nmax,atom->dihedral_per_atom); + if (atom->memcheck("dihedral_atom2")) + bytes += memory->usage(dihedral_atom2,nmax,atom->dihedral_per_atom); + if (atom->memcheck("dihedral_atom3")) + bytes += memory->usage(dihedral_atom3,nmax,atom->dihedral_per_atom); + if (atom->memcheck("dihedral_atom4")) + bytes += memory->usage(dihedral_atom4,nmax,atom->dihedral_per_atom); + + if (atom->memcheck("num_improper")) bytes += memory->usage(num_improper,nmax); + if (atom->memcheck("improper_type")) + bytes += memory->usage(improper_type,nmax,atom->improper_per_atom); + if (atom->memcheck("improper_atom1")) + bytes += memory->usage(improper_atom1,nmax,atom->improper_per_atom); + if (atom->memcheck("improper_atom2")) + bytes += memory->usage(improper_atom2,nmax,atom->improper_per_atom); + if (atom->memcheck("improper_atom3")) + bytes += memory->usage(improper_atom3,nmax,atom->improper_per_atom); + if (atom->memcheck("improper_atom4")) + bytes += memory->usage(improper_atom4,nmax,atom->improper_per_atom); + + return bytes; +} diff --git a/src/USER-GLE/atom_vec_gle.h b/src/USER-GLE/atom_vec_gle.h new file mode 100644 index 0000000000..d203cd4d2e --- /dev/null +++ b/src/USER-GLE/atom_vec_gle.h @@ -0,0 +1,99 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + * A class for an atom kind which carries around a variable number of + * additional degrees of freedom, to integrate a Generalized Langevin + * Equation dynamics. + * + * TODO: Probably this is better implemented as a hybrid type. + * ---------------------------------------------------------------------- */ + +#ifdef ATOM_CLASS + +AtomStyle(gle,AtomVecGLE) + +#else + +#ifndef LMP_ATOM_VEC_GLE_H +#define LMP_ATOM_VEC_GLE_H + +#include "atom_vec.h" + +namespace LAMMPS_NS { + +class AtomVecGLE : public AtomVec { + public: + AtomVecGLE(class LAMMPS *); //, int, char **); AtomVec classes used to accept arguments + virtual ~AtomVecGLE() {} + void grow(int); + void grow_reset(); + void copy(int, int, int); + virtual int pack_comm(int, int *, double *, int, int *); + virtual int pack_comm_vel(int, int *, double *, int, int *); + virtual void unpack_comm(int, int, double *); + virtual void unpack_comm_vel(int, int, double *); + int pack_reverse(int, int, double *); + void unpack_reverse(int, int *, double *); + virtual int pack_border(int, int *, double *, int, int *); + virtual int pack_border_vel(int, int *, double *, int, int *); + int pack_border_hybrid(int, int *, double *); + virtual void unpack_border(int, int, double *); + virtual void unpack_border_vel(int, int, double *); + int unpack_border_hybrid(int, int, double *); + virtual int pack_exchange(int, double *); + virtual int unpack_exchange(double *); + int size_restart(); + int pack_restart(int, double *); + int unpack_restart(double *); + void create_atom(int, double *); + void data_atom(double *, imageint, char **); + int data_atom_hybrid(int, char **); + void pack_data(double **); + int pack_data_hybrid(int, double *); + void write_data(FILE *, int, double **); + int write_data_hybrid(FILE *, double *); + bigint memory_usage(); + + protected: + int ns; + tagint *tag; + int *type,*mask; + imageint *image; + + // s[i][j] i is the atom index, j contains 3*(ns+1) items + // corresponding to x,y,z coordinates of different + // additional DOFs, i.e. s[0] will be [ x0 y0 z0 x1 y1 z1 ... ] + double **x, **v, **f, **s; + double *q; + tagint *molecule; + int **nspecial; + tagint **special; + int *num_bond; + int **bond_type; + tagint **bond_atom; + int *num_angle; + int **angle_type; + tagint **angle_atom1,**angle_atom2,**angle_atom3; + int *num_dihedral; + int **dihedral_type; + tagint **dihedral_atom1,**dihedral_atom2,**dihedral_atom3,**dihedral_atom4; + int *num_improper; + int **improper_type; + tagint **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4; +}; + +} + +#endif +#endif diff --git a/src/USER-GLE/fix_gle.cpp b/src/USER-GLE/fix_gle.cpp new file mode 100644 index 0000000000..91b55662dc --- /dev/null +++ b/src/USER-GLE/fix_gle.cpp @@ -0,0 +1,522 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Michele Ceriotti (Oxford), Joe Morrone (Columbia) +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "math.h" +#include "string.h" +#include "stdlib.h" +#include "fix_gle.h" +#include "math_extra.h" +#include "atom.h" +#include "atom_vec_ellipsoid.h" +#include "force.h" +#include "update.h" +#include "modify.h" +#include "compute.h" +#include "domain.h" +#include "region.h" +#include "respa.h" +#include "comm.h" +#include "input.h" +#include "variable.h" +#include "random_mars.h" +#include "memory.h" +#include "error.h" +#include "group.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{NOBIAS,BIAS}; +enum{CONSTANT,EQUAL,ATOM}; +//#define GLE_DEBUG 1 + +/* syntax for fix_gle: + * fix_gle nfix id-group gle ns amatrix [temp|cmatrix] seed + * + * */ + +/* ---------------------------------------------------------------------- */ + +namespace GLE { + /*!!MC!! here begins the GLE stuff*/ +#define midx(n,i,j) ((i)*(n)+(j)) + +//"stabilized" cholesky decomposition. does a LDL^t decomposition, then sets to zero the negative diagonal elements and gets MM^t +void StabCholesky(unsigned long n, const double* MMt, double* M) +{ + double L[n*n], D[n]; + + unsigned long i,j,k; + for(i=0; i0.?sqrt(D[i]):0.); + + for(i=0; i=0; p--) + { + MyMult(n, n, n, SM, EM, TMP); for (int i=0; iall(FLERR,"Illegal fix gle command"); + + ns = atof(arg[3]); + + if (ns>atom->ns) error->all(FLERR, "Atom kind has not enough GLE slots"); + + + A = new double[(ns+1)*(ns+1)]; + C = new double[(ns+1)*(ns+1)]; + T=S=gle_tmp=gle_rnd=NULL; gle_buff=0; + + // LOADS A matrix + printf("Reading A from %s\n", arg[4]); + FILE* fgle = fopen(arg[4], "r"); + + if (fgle==NULL) error->all(FLERR, "Cannot open A matrix GLE input"); + for (int i=0; i<(ns+1)*(ns+1); ++i) //ASSUMES A IS IN THE CORRECT INVERSE TIME UNITS + if (!fscanf(fgle,"%lf",&(A[i]))) error->all(FLERR, "Cannot read A matrix GLE input"); + + fclose(fgle); + + temp=atof(arg[5]); + if (temp==0.0) + { + printf("Reading C from %s\n", arg[5]); + fgle = fopen(arg[5], "r"); + if (fgle==NULL) error->all(FLERR, "Cannot open C matrix GLE input"); + for (int i=0; i<(ns+1)*(ns+1); ++i) //ASSUMES C IS IN THE CORRECT TEMPERATURE UNITS + if (!fscanf(fgle,"%lf",&(C[i]))) error->all(FLERR, "Cannot read C matrix GLE input"); + } + else + { + for (int i=0; i<(ns+1)*(ns+1); ++i) C[i]=0.0; + for (int i=0; i<(ns+1)*(ns+1); i+=(ns+2)) C[i]=temp*force->boltz/force->mvv2e; + } + +#ifdef GLE_DEBUG + printf("A Matrix\n"); + GLE::MyPrint(ns+1,A); + printf("C Matrix\n"); + GLE::MyPrint(ns+1,C); +#endif + + tau=10; + int seed = atoi(arg[6]); + + // initialize Marsaglia RNG with processor-unique seed + // NB: this means runs will not be the same with different numbers of processors + random = new RanMars(lmp,seed + comm->me); + + if (seed <= 0) error->all(FLERR,"Illegal fix gle command"); + + // allocate per-type arrays for mass-scaling + sqrt_m = new double[atom->ntypes+1]; + + energy = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +FixGLE::~FixGLE() +{ + delete random; + delete [] sqrt_m; + delete [] A; delete [] C; delete [] S; delete [] T; + delete [] gle_tmp; delete [] gle_rnd; +} + +/* ---------------------------------------------------------------------- */ + +int FixGLE::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + mask |= INITIAL_INTEGRATE_RESPA; + mask |= FINAL_INTEGRATE_RESPA; + mask |= THERMO_ENERGY; + + + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixGLE::init() +{ + + dogle = 1; + dorattle = 1; + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + + // compute Langevin terms + T = new double[(ns+1)*(ns+1)]; S = new double[(ns+1)*(ns+1)]; + double *tmp1 = new double[(ns+1)*(ns+1)], *tmp2 = new double[(ns+1)*(ns+1)]; + for (int i=0; i<(ns+1)*(ns+1); ++i) { tmp1[i]=-A[i]*update->dt*0.5; tmp2[i]=S[i]=0.0; } + GLE::MatrixExp(ns+1,tmp1,T); + + GLE::MyMult(ns+1,ns+1,ns+1,T,C,tmp1); + GLE::MyTrans(ns+1,T,tmp2); + GLE::MyMult(ns+1,ns+1,ns+1,tmp1,tmp2,S); + for (int i=0; i<(ns+1)*(ns+1); ++i) tmp1[i]=C[i]-S[i]; + + GLE::StabCholesky(ns+1, tmp1, S); //!TODO use symmetric square root, which is more stable. + +#ifdef GLE_DEBUG + printf("T Matrix\n"); + GLE::MyPrint(ns+1,T); + printf("S Matrix\n"); + GLE::MyPrint(ns+1,S); +#endif + + delete[] tmp1; delete[] tmp2; + + // set force prefactors + if (!atom->rmass) { + for (int i = 1; i <= atom->ntypes; i++) { + sqrt_m[i] = sqrt(atom->mass[i]); + } + } + + if (strstr(update->integrate_style,"respa")) { + nlevels_respa = ((Respa *) update->integrate)->nlevels; + step_respa = ((Respa *) update->integrate)->step; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixGLE::setup(int vflag) +{ + if (strstr(update->integrate_style,"verlet")) + post_force(vflag); + else { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(vflag,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + } +} + +void FixGLE::gle_integrate() +{ + double **v = atom->v; + double **s = atom->s; + double *rmass = atom->rmass, smi, ismi; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + + + if (nlocal > gle_buff) + { + if (gle_tmp !=NULL) { + delete [] gle_tmp; delete [] gle_rnd; + } + gle_tmp = new double[nlocal*(ns+1)]; + gle_rnd = new double[nlocal*(ns+1)]; + gle_buff = nlocal; + } + +#ifdef GLE_DEBUG + printf("!MC! GLE THERMO STEP dt=%f\n",update->dt); +#endif + + // gets kinetic energy before doing anything to the velocities + double deltae=0.0; + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ismi=rmass[i]; + for (int j = 0; j<3; ++j) { deltae-=ismi*v[i][j]*v[i][j]; } + } } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ismi=atom->mass[type[i]]; + for (int j = 0; j<3; ++j) { deltae-=ismi*v[i][j]*v[i][j]; } + } } + } + + // s is just taken to be the mass-scaled velocity! + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + smi=sqrt(rmass[i]); + for (int j = 0; j<3; ++j) { s[i][j] = v[i][j]*smi; } + } } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + smi=sqrt_m[type[i]]; + for (int j = 0; j<3; ++j) { s[i][j] = v[i][j]*smi; } + } } + } + + + for (int k = 0; k<3; k++) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) for (int j = 0; j<(ns+1); ++j) gle_rnd[j*(nlocal)+i] = s[i][3*j+k]; + } + GLE::MyMult(ns+1,nlocal,ns+1,T,gle_rnd,gle_tmp); + for (int i = 0; i < nlocal*(ns+1); i++) gle_rnd[i] = random->gaussian(); + GLE::MyMult(ns+1,nlocal,ns+1,S,gle_rnd,gle_tmp,1.0); + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) for (int j = 0; j<(ns+1); ++j) s[i][3*j+k] = gle_tmp[j*(nlocal)+i]; + } + } + + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ismi=1.0/sqrt(rmass[i]); + for (int j = 0; j<3; ++j) { v[i][j] = s[i][j]*ismi; } + } } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ismi=1.0/sqrt_m[type[i]]; + for (int j = 0; j<3; ++j) { v[i][j] = s[i][j]*ismi; } + } } + } + + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ismi=rmass[i]; + for (int j = 0; j<3; ++j) { deltae+=ismi*v[i][j]*v[i][j]; } + } } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ismi=atom->mass[type[i]]; + for (int j = 0; j<3; ++j) { deltae+=ismi*v[i][j]*v[i][j]; } + } } + } + + energy += deltae*0.5*force->mvv2e; +} + +void FixGLE::initial_integrate(int vflag) +{ + double dtfm; + + // update v and x of atoms in group + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + + if (dogle) gle_integrate(); + +#ifdef GLE_DEBUG + printf("!MC! GLE P1 STEP dt=%f\n",dtv); + printf("!MC! GLE Q STEP\n"); +#endif + + if (rmass) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } +} + +void FixGLE::final_integrate() +{ + double dtfm; + + // update v of atoms in group + + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + if (rmass) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + } + + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + } + } + +#ifdef GLE_DEBUG + printf("!MC! GLE P2 STEP dt=%f\n",dtv); +#endif + if (dogle) gle_integrate(); +} +/* ---------------------------------------------------------------------- */ + +void FixGLE::initial_integrate_respa(int vflag, int ilevel, int iloop) +{ + dtv = step_respa[ilevel]; + dtf = 0.5 * step_respa[ilevel] * force->ftm2v; + + // innermost level - NVE update of v and x + // all other levels - NVE update of v + + if (ilevel==nlevels_respa-1) gle_integrate(); + dogle=0; + if (ilevel == 0) initial_integrate(vflag); + else { dorattle=1; final_integrate();} // If RATTLE S2 should not be called here, dorattle=0 +} + +void FixGLE::final_integrate_respa(int ilevel, int iloop) +{ + dtv = step_respa[ilevel]; //!MC! + dtf = 0.5 * step_respa[ilevel] * force->ftm2v; + dogle=0; + if(0) { // If true turn on stage 2 rattle for final integrate only + if(ilevel==nlevels_respa-1) + dorattle=1; + else + dorattle=0; + } else { dorattle=1; } + final_integrate(); + if (ilevel==nlevels_respa-1) gle_integrate(); +} + + +double FixGLE::compute_scalar() +{ + double energy_me = energy; + double energy_all; + MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); + + return -energy_all; +} + +/* ---------------------------------------------------------------------- + extract thermostat properties +------------------------------------------------------------------------- */ + +void *FixGLE::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"t_target") == 0) { + return &temp; + } + return NULL; +} + diff --git a/src/USER-GLE/fix_gle.h b/src/USER-GLE/fix_gle.h new file mode 100644 index 0000000000..b715d86d1c --- /dev/null +++ b/src/USER-GLE/fix_gle.h @@ -0,0 +1,61 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(gle,FixGLE) + +#else + +#ifndef LMP_FIX_GLE_H +#define LMP_FIX_GLE_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixGLE : public Fix { + public: + FixGLE(class LAMMPS *, int, char **); + virtual ~FixGLE(); + int setmask(); + void init(); + void setup(int); + void gle_integrate(); + void initial_integrate_respa(int vflag, int ilevel, int iloop); + void final_integrate_respa(int ilevel, int iloop); + void initial_integrate(int vflag); + void final_integrate(); + double compute_scalar(); + virtual void *extract(const char *, int &); + + protected: + int ns; + double *A, *C, *S, *T; + double *gle_tmp, *gle_rnd; int gle_buff; + double temp, tau, dtv, dtf; + int dogle; + int dorattle; + class RanMars *random; + double *sqrt_m; + double *step_respa; + double energy; + int nlevels_respa; + + double **vaux; +}; + +} + +#endif +#endif diff --git a/src/atom.cpp b/src/atom.cpp index 5554117f5c..ff2b15f4a0 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -52,6 +52,7 @@ using namespace MathConst; Atom::Atom(LAMMPS *lmp) : Pointers(lmp) { natoms = 0; + ns = 0; //!GLE nlocal = nghost = nmax = 0; ntypes = 0; nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0; @@ -72,6 +73,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) type = mask = NULL; image = NULL; x = v = f = NULL; + s = NULL; //!GLE molecule = NULL; molindex = molatom = NULL; @@ -205,6 +207,7 @@ Atom::~Atom() memory->destroy(x); memory->destroy(v); memory->destroy(f); + memory->destroy(s); //!GLE memory->destroy(molecule); memory->destroy(molindex); @@ -1891,6 +1894,7 @@ void *Atom::extract(char *name) if (strcmp(name,"x") == 0) return (void *) x; if (strcmp(name,"v") == 0) return (void *) v; if (strcmp(name,"f") == 0) return (void *) f; + if (strcmp(name,"s") == 0) return (void *) s; //!GLE if (strcmp(name,"molecule") == 0) return (void *) molecule; if (strcmp(name,"q") == 0) return (void *) q; if (strcmp(name,"mu") == 0) return (void *) mu; diff --git a/src/atom.h b/src/atom.h index c6bebe88a9..882563cf55 100644 --- a/src/atom.h +++ b/src/atom.h @@ -51,6 +51,7 @@ class Atom : protected Pointers { imageint *image; double **x,**v,**f; + double **s; int ns; //!GLE tagint *molecule; int *molindex,*molatom;