Langevin Simulations in 1D

Hi,

I’m trying to simulate a 1D model of a “wave” where there’s only motion in the y dimension and there are a lot of fix particles. Nevertheless, when I’m running a Langevin thermostat (Which should only work in 1D), I’m getting a much lower temperature than I should. Do you notice what I’m doing wrong and how to fix it?

Test system for LAMMPS simulation

units real
atom_style full
log sim438.log

For a single processor calculation

variable T equal 273 # Simulation temperature
variable salt equal 100.0 # Salt concentration [mM]

Random number seed for Langevin integrator

variable random equal 12345

Specify the different interaction styles

bond_style hybrid morse force2
angle_style none
dihedral_style none
pair_style none

Periodic boundary conditions

boundary p p p

Turn on Newton’s 2nd law

newton on #yes

Read in the configuration

read_data bdna_conf.in

Specify parameters for the neighbor list

neighbor 4.0 multi
neigh_modify check no
atom_modify sort 0 0.0

A timestep of 0.01 ps

timestep 10

Define the groups of the real atoms and their fix initial positions

group realatom1 id 2:10 24:32 46:54 68:76 90:98 112:120 134:142 156:164 178:186 200:208 222:230 244:252 266:274 288:296 310:318 332:340 354:362 376:384 398:406 420:428 442:450 464:472 486:494 508:516 530:538 552:560 574:582 596:604 618:626 640:648 662:670 684:692 706:714 728:736 750:758 772:780 794:802 816:824 838:846 860:868 882:890 904:912 926:934 948:956 970:978 992:1000 1014:1022 1036:1044 1058:1066 1080:1088 1102:1110 1124:1132 1146:1154 1168:1176 1190:1198 1212:1220 1234:1242 1256:1264 1278:1286 1300:1308 1322:1330 1344:1352 1366:1374 1388:1396 1410:1418 1432:1440 1454:1462 1476:1484 1498:1506 1520:1528 1542:1550 1564:1572 1586:1594 1608:1616 1630:1638 1652:1660 1674:1682 1696:1704 1718:1726 1740:1748 1762:1770 1784:1792 1806:1814 1828:1836 1850:1858 1872:1880 1894:1902 1916:1924 1938:1946 1960:1968 1982:1990 2004:2012 2026:2034 2048:2056 2070:2078 2092:2100 2114:2122 2136:2144 2158:2166 2180:2188 2202:2210 2224:2232 2246:2254 2268:2276 2290:2298 2312:2320 2334:2342 2356:2364 2378:2386 2400:2408 2422:2430 2444:2452 2466:2474 2488:2496 2510:2518 2532:2540 2554:2562 2576:2584 2598:2606 2620:2628 2642:2650 2664:2672 2686:2694 2708:2716 2730:2738 2752:2760 2774:2782 2796:2804 2818:2826 2840:2848 2862:2870 2884:2892 2906:2914 2928:2936 2950:2958 2972:2980 2994:3002 3016:3024 3038:3046 3060:3068 3082:3090 3104:3112 3126:3134 3148:3156 3170:3178 3192:3200 3214:3222 3236:3244 3258:3266 3280:3288 3302:3310 3324:3332 3346:3354 3368:3376 3390:3398 3412:3420 3434:3442 3456:3464 3478:3486 3500:3508 3522:3530
group fixatom1 type 5
group boundaries id 1 11 23 33 45 55 67 77 89 99 111 121 133 143 155 165 177 187 199 209 221 231 243 253 265 275 287 297 309 319 331 341 353 363 375 385 397 407 419 429 441 451 463 473 485 495 507 517 529 539 551 561 573 583 595 605 617 627 639 649 661 671 683 693 705 715 727 737 749 759 771 781 793 803 815 825 837 847 859 869 881 891 903 913 925 935 947 957 969 979 991 1001 1013 1023 1035 1045 1057 1067 1079 1089 1101 1111 1123 1133 1145 1155 1167 1177 1189 1199 1211 1221 1233 1243 1255 1265 1277 1287 1299 1309 1321 1331 1343 1353 1365 1375 1387 1397 1409 1419 1431 1441 1453 1463 1475 1485 1497 1507 1519 1529 1541 1551 1563 1573 1585 1595 1607 1617 1629 1639 1651 1661 1673 1683 1695 1705 1717 1727 1739 1749 1761 1771 1783 1793 1805 1815 1827 1837 1849 1859 1871 1881 1893 1903 1915 1925 1937 1947 1959 1969 1981 1991 2003 2013 2025 2035 2047 2057 2069 2079 2091 2101 2113 2123 2135 2145 2157 2167 2179 2189 2201 2211 2223 2233 2245 2255 2267 2277 2289 2299 2311 2321 2333 2343 2355 2365 2377 2387 2399 2409 2421 2431 2443 2453 2465 2475 2487 2497 2509 2519 2531 2541 2553 2563 2575 2585 2597 2607 2619 2629 2641 2651 2663 2673 2685 2695 2707 2717 2729 2739 2751 2761 2773 2783 2795 2805 2817 2827 2839 2849 2861 2871 2883 2893 2905 2915 2927 2937 2949 2959 2971 2981 2993 3003 3015 3025 3037 3047 3059 3069 3081 3091 3103 3113 3125 3135 3147 3157 3169 3179 3191 3201 3213 3223 3235 3245 3257 3267 3279 3289 3301 3311 3323 3333 3345 3355 3367 3377 3389 3399 3411 3421 3433 3443 3455 3465 3477 3487 3499 3509 3521 3531

Initialize velocities from a Gaussian distribution

velocity realatom1 create {T} {random} rot yes mom yes dist gaussian
velocity realatom1 set 0.0 NULL 0.0
velocity fixatom1 set 0.0 0.0 0.0
velocity boundaries set 0.0 0.0 0.0
run 0
velocity all scale ${T}

Calculating the different components of the non-bonded energy

compute 1 all bond/local dist eng
compute newtemp all temp/partial 0 1 0

Specifying the frequency of thermodynamic output

thermo 100
thermo_style custom step ebond temp

Specifying a Langevin integrator to perform a simulation in the NVT ensemble

fix 1 realatom1 langevin {T} {T} 5000 ${random} gjf yes
fix 2 realatom1 nve
fix 3 realatom1 setforce 0.0 NULL 0.0
fix 4 fixatom1 setforce 0.0 0.0 0.0
fix 5 boundaries setforce 0.0 0.0 0.0
fix_modify 1 temp newtemp

Write configuration to file

dump 1 realatom1 custom 100 traj438.xyz x y z vy

Run X number of steps

run 30000

You have fixed (immobilized) atoms/particles, so how are you computing and outputting the temperature? Is it a temperature of the entire system or just of the free (mobile) atoms?

Ray

Dear Ray,

thanks for your reply! I’m computing the temperature of only the mobile ones.

I don’t think so…

This command "compute newtemp all temp/partial 0 1 0 " is applied to group all, and the thermo_style outputs “temp”, which is also group all with all three components included.

Ray

Hi Ray, thanks a lot for your reply!!! I tried changing it to what follows but I still get a lower temperature than I should,

what would you change?

I tried to write “thermo_style custom step ebond newtemp” but it gives me an error

Calculating the different components of the non-bonded energy

compute 1 all bond/local dist eng
compute newtemp realatom1 temp/partial 0 1 0

Specifying the frequency of thermodynamic output

thermo 500
thermo_style custom step ebond temp

Specifying a Langevin integrator to perform a simulation in the NVT ensemble

fix 3 realatom1 langevin {T} {T} 5000 ${random} gjf yes
fix 4 realatom1 nve
fix 5 realatom1 setforce 0.0 NULL 0.0
fix 6 fixatom1 setforce 0.0 0.0 0.0
fix 7 boundaries setforce 0.0 0.0 0.0
fix_modify 3 temp newtemp

Write configuration to file

dump 1 mobile custom 500 traj273.xyz x y z vy

Hi Ray, thanks a lot for your reply!!! I tried changing it to what follows
but I still get a lower temperature than I should,

Of course you did since thermo_style is still outputting "temp" ...

what would you change?
I tried to write "thermo_style custom step ebond newtemp" but it gives me
an error

Since it is from a compute you will have to use "c_newtemp". See (
LAMMPS Molecular Dynamics Simulator) for more details.

Ray

You were right, Ray! It worked!!! Thank you very much!

There still seems to be something odd with the thermostat. This is my current input, and I’m getting way lower displacements than I expected…

Test system for LAMMPS simulation

units real
atom_style full
log sim323.log

For a single processor calculation

variable T equal 423 # Simulation temperature
variable salt equal 100.0 # Salt concentration [mM]

Random number seed for Langevin integrator

variable random equal 12345

Specify the different interaction styles

bond_style hybrid morse potential1
angle_style none
dihedral_style none
pair_style none

Periodic boundary conditions

boundary p p p

Turn on Newton’s 2nd law

newton on #yes

Read in the configuration

read_data bdna_conf.in

Specify parameters for the neighbor list

neighbor 4.0 multi
neigh_modify check no
atom_modify sort 0 0.0

A timestep of 0.10 ps

timestep 100

Define the groups of the real atoms and their fix initial positions

group realatom1 id 2:10 24:32 46:54 68:76 90:98 112:120 134:142 156:164 178:186 200:208 222:230 244:252 266:274 288:296 310:318 332:340 354:362 376:384 398:406 420:428 442:450 464:472 486:494 508:516 530:538 552:560 574:582 596:604 618:626 640:648 662:670 684:692 706:714 728:736 750:758 772:780 794:802 816:824 838:846 860:868 882:890 904:912 926:934 948:956 970:978 992:1000 1014:1022 1036:1044 1058:1066 1080:1088 1102:1110 1124:1132 1146:1154 1168:1176 1190:1198 1212:1220 1234:1242 1256:1264 1278:1286 1300:1308 1322:1330 1344:1352 1366:1374 1388:1396 1410:1418 1432:1440 1454:1462 1476:1484 1498:1506 1520:1528 1542:1550 1564:1572 1586:1594 1608:1616 1630:1638 1652:1660 1674:1682 1696:1704 1718:1726 1740:1748 1762:1770 1784:1792 1806:1814 1828:1836 1850:1858 1872:1880 1894:1902 1916:1924 1938:1946 1960:1968 1982:1990 2004:2012 2026:2034 2048:2056 2070:2078 2092:2100 2114:2122 2136:2144 2158:2166 2180:2188 2202:2210 2224:2232 2246:2254 2268:2276 2290:2298 2312:2320 2334:2342 2356:2364 2378:2386 2400:2408 2422:2430 2444:2452 2466:2474 2488:2496 2510:2518 2532:2540 2554:2562 2576:2584 2598:2606 2620:2628 2642:2650 2664:2672 2686:2694 2708:2716 2730:2738 2752:2760 2774:2782 2796:2804 2818:2826 2840:2848 2862:2870 2884:2892 2906:2914 2928:2936 2950:2958 2972:2980 2994:3002 3016:3024 3038:3046 3060:3068 3082:3090 3104:3112 3126:3134 3148:3156 3170:3178 3192:3200 3214:3222 3236:3244 3258:3266 3280:3288 3302:3310 3324:3332 3346:3354 3368:3376 3390:3398 3412:3420 3434:3442 3456:3464 3478:3486 3500:3508 3522:3530
group fixatom1 type 5
group boundaries id 1 11 23 33 45 55 67 77 89 99 111 121 133 143 155 165 177 187 199 209 221 231 243 253 265 275 287 297 309 319 331 341 353 363 375 385 397 407 419 429 441 451 463 473 485 495 507 517 529 539 551 561 573 583 595 605 617 627 639 649 661 671 683 693 705 715 727 737 749 759 771 781 793 803 815 825 837 847 859 869 881 891 903 913 925 935 947 957 969 979 991 1001 1013 1023 1035 1045 1057 1067 1079 1089 1101 1111 1123 1133 1145 1155 1167 1177 1189 1199 1211 1221 1233 1243 1255 1265 1277 1287 1299 1309 1321 1331 1343 1353 1365 1375 1387 1397 1409 1419 1431 1441 1453 1463 1475 1485 1497 1507 1519 1529 1541 1551 1563 1573 1585 1595 1607 1617 1629 1639 1651 1661 1673 1683 1695 1705 1717 1727 1739 1749 1761 1771 1783 1793 1805 1815 1827 1837 1849 1859 1871 1881 1893 1903 1915 1925 1937 1947 1959 1969 1981 1991 2003 2013 2025 2035 2047 2057 2069 2079 2091 2101 2113 2123 2135 2145 2157 2167 2179 2189 2201 2211 2223 2233 2245 2255 2267 2277 2289 2299 2311 2321 2333 2343 2355 2365 2377 2387 2399 2409 2421 2431 2443 2453 2465 2475 2487 2497 2509 2519 2531 2541 2553 2563 2575 2585 2597 2607 2619 2629 2641 2651 2663 2673 2685 2695 2707 2717 2729 2739 2751 2761 2773 2783 2795 2805 2817 2827 2839 2849 2861 2871 2883 2893 2905 2915 2927 2937 2949 2959 2971 2981 2993 3003 3015 3025 3037 3047 3059 3069 3081 3091 3103 3113 3125 3135 3147 3157 3169 3179 3191 3201 3213 3223 3235 3245 3257 3267 3279 3289 3301 3311 3323 3333 3345 3355 3367 3377 3389 3399 3411 3421 3433 3443 3455 3465 3477 3487 3499 3509 3521 3531
group mobile union realatom1 boundaries

Initialize velocities from a Gaussian distribution

velocity mobile create {T} {random} rot yes mom yes dist gaussian
velocity mobile set 0.0 NULL 0.0
velocity fixatom1 set 0.0 0.0 0.0
run 0
velocity all scale ${T}

Calculating the different components of the non-bonded energy

compute 1 all bond/local dist eng
compute newtemp mobile temp/partial 0 1 0

Specifying the frequency of thermodynamic output

thermo 500
thermo_style custom step ebond etotal c_newtemp

Specifying a Langevin integrator to perform a simulation in the NVT ensemble

fix 3 mobile langevin {T} {T} 5000 ${random} gjf yes
fix 4 mobile nve
fix 5 mobile setforce 0.0 NULL 0.0
fix 6 fixatom1 setforce 0.0 0.0 0.0

Write configuration to file

dump 1 mobile custom 500 traj323.xyz x y z vy
dump 2 all xyz 500 traj323b.xyz

Run X number of steps

run 30000

Well, “my script does not do what I wanted it to do” is not exactly how this mailing list works…

A detailed description of the problem you observed, followed by pasting only the relevant commands (but a complete script that runs) is the best way to go.

Ray

You’re right. I think there is a problem with the thermostat, which is either

(1) working at 3D and then making the other coordinates zero, which will lead to low displacements, or

(2) working at a lower temperature than it should, maybe due to the fixed atoms

How did you analyze/estimate the displacement and came to the conclusion that “displacements” are lower(smaller)? And what exactly did you mean by “displacements”? Is it a averaged or per-atom property?

These questions are becoming more about your model and its analysis, which is drifting away from the purpose of this mailing list. I suggest you discuss your hypotheses with your advisor and/or colleagues as it is also not the purpose here. Hope you understand.

Ray

Thanks! I’m running the simulations at a lot of temperatures and seeing how much each atom moves and they just don’t move much when working with very high temperatures, studying at a per-atom level.

I have coded a 1D langevin thermostat previously and it is providing higher movements in comparison, that’s why I think there has to be something wrong with this method.

I could work straight with the 1D thermostat, however, I’d like to learn if there is a way of using the usual langevin thermostat and making it act in only 1D, as I haven’t seen anything like that discussed in the mailing list and seems quite convenient (in 2D there is the fix enforce2d command to help with this kind of issues)