@@ -59,10 +59,11 @@ def write_lammpsdata(
59
59
'full', 'atomic', 'charge', 'molecular'
60
60
see http://lammps.sandia.gov/doc/atom_style.html for more information
61
61
on atom styles.
62
- unit_style : str, optional, default='real'
63
- Defines to unit style to be save in a LAMMPS data file. Current styles
64
- are supported: 'real', 'lj', see lammps unit style documentation:
65
- https://lammps.sandia.gov/doc/99/units.html for more information.
62
+ unit_style: str, optional, default='real'
63
+ Defines to unit style to be save in a LAMMPS data file. Defaults to
64
+ 'real' units. Current styles are supported: 'real', 'lj', 'metal'
65
+ see https://lammps.sandia.gov/doc/99/units.html for more information
66
+ on unit styles
66
67
mins : list, optional, default=None
67
68
Minimum box dimension in x, y, z directions, nm
68
69
maxs : list, optional, default=None
@@ -162,7 +163,7 @@ def write_lammpsdata(
162
163
atom_style
163
164
)
164
165
)
165
- if unit_style not in ["real" , "lj" ]:
166
+ if unit_style not in ["real" , "lj" , "metal" ]:
166
167
raise ValueError (
167
168
'Unit style "{}" is invalid or is not currently supported' .format (
168
169
unit_style
@@ -305,10 +306,23 @@ def write_lammpsdata(
305
306
* 10 ** - 6
306
307
)
307
308
charges [np .isinf (charges )] = 0
308
- else :
309
+ eng_unit_str = " "
310
+
311
+ elif unit_style == "real" :
309
312
sigma_conversion_factor = 1
310
313
epsilon_conversion_factor = 1
311
314
mass_conversion_factor = 1
315
+ eng_unit_str = "kcal/mol"
316
+
317
+ elif unit_style == "metal" :
318
+ sigma_conversion_factor = 1 # unit in Angstrom
319
+ epsilon_conversion_factor = 1 / 0.043364115 # kcal/mol to eV
320
+ mass_conversion_factor = 1
321
+ eng_unit_str = "eV"
322
+ else :
323
+ raise ValueError (
324
+ "Currently only support 'real', 'lj', and 'metal' unit_styles."
325
+ )
312
326
313
327
# Divide by conversion factor
314
328
Lx = box .Lx * (1 / sigma_conversion_factor )
@@ -441,12 +455,13 @@ def write_lammpsdata(
441
455
pair_coeff_label ,
442
456
unit_style ,
443
457
nbfix_in_data_file ,
458
+ eng_unit_str ,
444
459
)
445
460
446
461
# Write bond coefficients
447
462
if bonds :
448
463
_write_bond_information (
449
- structure , data , unique_bond_types , unit_style
464
+ structure , data , unique_bond_types , unit_style , eng_unit_str
450
465
)
451
466
452
467
# Write angle coefficients
@@ -457,6 +472,7 @@ def write_lammpsdata(
457
472
unique_angle_types ,
458
473
use_urey_bradleys ,
459
474
unit_style ,
475
+ eng_unit_str ,
460
476
)
461
477
462
478
# Write dihedral coefficients
@@ -468,6 +484,7 @@ def write_lammpsdata(
468
484
unit_style ,
469
485
use_rb_torsions ,
470
486
use_dihedrals ,
487
+ eng_unit_str ,
471
488
)
472
489
473
490
# Write improper coefficients
@@ -487,14 +504,14 @@ def write_lammpsdata(
487
504
if atom_style == "atomic" :
488
505
atom_line = "{index:d}\t {type_index:d}\t {x:.6f}\t {y:.6f}\t {z:.6f}\n "
489
506
elif atom_style == "charge" :
490
- if unit_style == "real" :
507
+ if unit_style in [ "real" , "metal" ] :
491
508
atom_line = "{index:d}\t {type_index:d}\t {charge:.6f}\t {x:.6f}\t {y:.6f}\t {z:.6f}\n "
492
509
elif unit_style == "lj" :
493
510
atom_line = "{index:d}\t {type_index:d}\t {charge:.4ef}\t {x:.6f}\t {y:.6f}\t {z:.6f}\n "
494
511
elif atom_style == "molecular" :
495
512
atom_line = "{index:d}\t {zero:d}\t {type_index:d}\t {x:.6f}\t {y:.6f}\t {z:.6f}\n "
496
513
elif atom_style == "full" :
497
- if unit_style == "real" :
514
+ if unit_style in [ "real" , "metal" ] :
498
515
atom_line = "{index:d}\t {zero:d}\t {type_index:d}\t {charge:.6f}\t {x:.6f}\t {y:.6f}\t {z:.6f}\n "
499
516
elif unit_style == "lj" :
500
517
atom_line = "{index:d}\t {zero:d}\t {type_index:d}\t {charge:.4e}\t {x:.6f}\t {y:.6f}\t {z:.6f}\n "
@@ -1167,6 +1184,7 @@ def _write_pair_information(
1167
1184
pair_coeff_label ,
1168
1185
unit_style ,
1169
1186
nbfix_in_data_file ,
1187
+ eng_unit_str ,
1170
1188
):
1171
1189
"""Write nonbonded pair information to lammps data file."""
1172
1190
epsilons = (
@@ -1255,7 +1273,9 @@ def _write_pair_information(
1255
1273
else :
1256
1274
data .write ("\n PairIJ Coeffs # modified lj\n " )
1257
1275
1258
- data .write ("# type1 type2\t epsilon (kcal/mol)\t sigma (Angstrom)\n " )
1276
+ data .write (
1277
+ "# type1 type2\t epsilon (%s)\t sigma (Angstrom)\n " % eng_unit_str
1278
+ )
1259
1279
1260
1280
for (type1 , type2 ), (sigma , epsilon ) in coeffs .items ():
1261
1281
data .write (
@@ -1279,7 +1299,9 @@ def _write_pair_information(
1279
1299
"{}\t {:.5f}\t {:.5f}\n " .format (idx , epsilon , sigma_dict [idx ])
1280
1300
)
1281
1301
print ("Copy these commands into your input script:\n " )
1282
- print ("# type1 type2\t epsilon (kcal/mol)\t sigma (Angstrom)\n " )
1302
+ print (
1303
+ "# type1 type2\t epsilon (%s)\t sigma (Angstrom)\n " % eng_unit_str
1304
+ )
1283
1305
for (type1 , type2 ), (sigma , epsilon ) in coeffs .items ():
1284
1306
print (
1285
1307
"pair_coeff\t {0} \t {1} \t {2} \t \t {3} \t \t # {4} \t {5}" .format (
@@ -1299,8 +1321,8 @@ def _write_pair_information(
1299
1321
else :
1300
1322
data .write ("\n Pair Coeffs # lj\n " )
1301
1323
1302
- if unit_style == "real" :
1303
- data .write ("#\t epsilon (kcal/mol )\t \t sigma (Angstrom)\n " )
1324
+ if unit_style in [ "real" , "metal" ] :
1325
+ data .write ("#\t epsilon (%s )\t \t sigma (Angstrom)\n " % eng_unit_str )
1304
1326
elif unit_style == "lj" :
1305
1327
data .write ("#\t reduced_epsilon \t \t reduced_sigma \n " )
1306
1328
for idx , epsilon in sorted (epsilon_dict .items ()):
@@ -1311,11 +1333,13 @@ def _write_pair_information(
1311
1333
)
1312
1334
1313
1335
1314
- def _write_bond_information (structure , data , unique_bond_types , unit_style ):
1336
+ def _write_bond_information (
1337
+ structure , data , unique_bond_types , unit_style , eng_unit_str
1338
+ ):
1315
1339
"""Write Bond Coeffs section of lammps data file."""
1316
1340
data .write ("\n Bond Coeffs # harmonic\n " )
1317
- if unit_style == "real" :
1318
- data .write ("#\t k(kcal/mol/ angstrom^2)\t \t req(angstrom)\n " )
1341
+ if unit_style == "real" or unit_style == "metal" :
1342
+ data .write ("#\t k(%s/ angstrom^2)\t \t req(angstrom)\n " % eng_unit_str )
1319
1343
elif unit_style == "lj" :
1320
1344
data .write ("#\t reduced_k\t \t reduced_req\n " )
1321
1345
sorted_bond_types = {
@@ -1335,7 +1359,12 @@ def _write_bond_information(structure, data, unique_bond_types, unit_style):
1335
1359
1336
1360
1337
1361
def _write_angle_information (
1338
- structure , data , unique_angle_types , use_urey_bradleys , unit_style
1362
+ structure ,
1363
+ data ,
1364
+ unique_angle_types ,
1365
+ use_urey_bradleys ,
1366
+ unit_style ,
1367
+ eng_unit_str ,
1339
1368
):
1340
1369
"""Write Angle Coeffs section of lammps data file."""
1341
1370
sorted_angle_types = {
@@ -1345,7 +1374,8 @@ def _write_angle_information(
1345
1374
if use_urey_bradleys :
1346
1375
data .write ("\n Angle Coeffs # charmm\n " )
1347
1376
data .write (
1348
- "#\t k(kcal/mol/rad^2)\t \t theteq(deg)\t k(kcal/mol/angstrom^2)\t req(angstrom)\n "
1377
+ "#\t k(%s/rad^2)\t \t theteq(deg)\t k(%s/angstrom^2)\t req(angstrom)\n "
1378
+ % (eng_unit_str , eng_unit_str )
1349
1379
)
1350
1380
for params , idx in sorted_angle_types .items ():
1351
1381
data .write ("{}\t {}\t {:.5f}\t {:.5f}\t {:.5f}\n " .format (idx , * params ))
@@ -1356,7 +1386,7 @@ def _write_angle_information(
1356
1386
if unit_style == "lj" :
1357
1387
data .write ("#\t reduced_k\t \t theteq(deg)\n " )
1358
1388
else :
1359
- data .write ("#\t k(kcal/mol/ rad^2)\t \t theteq(deg)\n " )
1389
+ data .write ("#\t k(%s/ rad^2)\t \t theteq(deg)\n " % eng_unit_str )
1360
1390
1361
1391
for params , idx in sorted_angle_types .items ():
1362
1392
data .write (
@@ -1378,6 +1408,7 @@ def _write_dihedral_information(
1378
1408
unit_style ,
1379
1409
use_rb_torsions ,
1380
1410
use_dihedrals ,
1411
+ eng_unit_str ,
1381
1412
):
1382
1413
"""Write Dihedral Coeffs section of lammps data file."""
1383
1414
sorted_dihedral_types = {
@@ -1388,9 +1419,10 @@ def _write_dihedral_information(
1388
1419
}
1389
1420
if use_rb_torsions :
1390
1421
data .write ("\n Dihedral Coeffs # opls\n " )
1391
- if unit_style == "real" :
1422
+ if unit_style == "real" or unit_style == "metal" :
1392
1423
data .write (
1393
- "#\t f1(kcal/mol)\t f2(kcal/mol)\t f3(kcal/mol)\t f4(kcal/mol)\n "
1424
+ "#\t f1(%s)\t f2(%s)\t f3(%s)\t f4(%s)\n "
1425
+ % (eng_unit_str , eng_unit_str , eng_unit_str , eng_unit_str )
1394
1426
)
1395
1427
elif unit_style == "lj" :
1396
1428
data .write ("#\t f1\t f2\t f3\t f4 (all lj reduced units)\n " )
@@ -1422,11 +1454,11 @@ def _write_dihedral_information(
1422
1454
data .write ("#k, n, phi, weight\n " )
1423
1455
for params , idx in sorted_dihedral_types .items ():
1424
1456
data .write (
1425
- "{}\t {:.5f}\t {:d}\t {:.2f }\t {:.5f}\t # {}\t {}\t {}\t {}\n " .format (
1457
+ "{}\t {:.5f}\t {:d}\t {:d }\t {:.5f}\t # {}\t {}\t {}\t {}\n " .format (
1426
1458
idx ,
1427
1459
params [0 ],
1428
1460
params [1 ],
1429
- params [2 ],
1461
+ int ( params [2 ]) ,
1430
1462
params [3 ],
1431
1463
params [6 ],
1432
1464
params [7 ],
0 commit comments