SciCodePile/SciCode-Domain-Code
DATA1: Domain-Specific Code Dataset Dataset Overview DATA1 is a large-scale domain-specific code dataset focusing on code samples from interdisciplinary fields such as biology, chemistry, materials science, and related areas. The dataset is collected and organized from GitHub repositories, covering 178 different domain topics with over 1.1 billion lines of code. Dataset Statistics Total Datasets: 178 CSV files Total Data Size: ~115 GB Total Lines… See the full description on the dataset page: https://huggingface.co/datasets/SciCodePile/SciCode-Domain-Code.
42.4k
1"keyword","repo_name","file_path","file_extension","file_size","line_count","content","language"
2"Multi-scale modeling","vishalsubbiah/Grain-Boundary-Energies-LAMMPS","Code and Scripts/Python and Lammps/FullStackAll/FullStack101010/Experiments/Fe/300/fractions.py",".py","22391","607","# Originally contributed by Sjoerd Mullender.3# Significantly modified by Jeffrey Yasskin <jyasskin at gmail.com>.4 5""""""Rational, infinite-precision, real numbers.""""""6 7from __future__ import division8from decimal import Decimal9import math10import numbers11import operator12import re13 14__all__ = ['Fraction', 'gcd']15 16Rational = numbers.Rational17 18 19def gcd(a, b):20 """"""Calculate the Greatest Common Divisor of a and b.21 22 Unless b==0, the result will have the same sign as b (so that when23 b is divided by it, the result comes out positive).24 """"""25 while b:26 a, b = b, a%b27 return a28 29 30_RATIONAL_FORMAT = re.compile(r""""""31 \A\s* # optional whitespace at the start, then32 (?P<sign>[-+]?) # an optional sign, then33 (?=\d|\.\d) # lookahead for digit or .digit34 (?P<num>\d*) # numerator (possibly empty)35 (?: # followed by36 (?:/(?P<denom>\d+))? # an optional denominator37 | # or38 (?:\.(?P<decimal>\d*))? # an optional fractional part39 (?:E(?P<exp>[-+]?\d+))? # and optional exponent40 )41 \s*\Z # and optional whitespace to finish42"""""", re.VERBOSE | re.IGNORECASE)43 44 45class Fraction(Rational):46 """"""This class implements rational numbers.47 48 In the two-argument form of the constructor, Fraction(8, 6) will49 produce a rational number equivalent to 4/3. Both arguments must50 be Rational. The numerator defaults to 0 and the denominator51 defaults to 1 so that Fraction(3) == 3 and Fraction() == 0.52 53 Fractions can also be constructed from:54 55 - numeric strings similar to those accepted by the56 float constructor (for example, '-2.3' or '1e10')57 58 - strings of the form '123/456'59 60 - float and Decimal instances61 62 - other Rational instances (including integers)63 64 """"""65 66 __slots__ = ('_numerator', '_denominator')67 68 # We're immutable, so use __new__ not __init__69 def __new__(cls, numerator=0, denominator=None):70 """"""Constructs a Fraction.71 72 Takes a string like '3/2' or '1.5', another Rational instance, a73 numerator/denominator pair, or a float.74 75 Examples76 --------77 78 >>> Fraction(10, -8)79 Fraction(-5, 4)80 >>> Fraction(Fraction(1, 7), 5)81 Fraction(1, 35)82 >>> Fraction(Fraction(1, 7), Fraction(2, 3))83 Fraction(3, 14)84 >>> Fraction('314')85 Fraction(314, 1)86 >>> Fraction('-35/4')87 Fraction(-35, 4)88 >>> Fraction('3.1415') # conversion from numeric string89 Fraction(6283, 2000)90 >>> Fraction('-47e-2') # string may include a decimal exponent91 Fraction(-47, 100)92 >>> Fraction(1.47) # direct construction from float (exact conversion)93 Fraction(6620291452234629, 4503599627370496)94 >>> Fraction(2.25)95 Fraction(9, 4)96 >>> Fraction(Decimal('1.47'))97 Fraction(147, 100)98 99 """"""100 self = super(Fraction, cls).__new__(cls)101 102 if denominator is None:103 if isinstance(numerator, Rational):104 self._numerator = numerator.numerator105 self._denominator = numerator.denominator106 return self107 108 elif isinstance(numerator, float):109 # Exact conversion from float110 value = Fraction.from_float(numerator)111 self._numerator = value._numerator112 self._denominator = value._denominator113 return self114 115 elif isinstance(numerator, Decimal):116 value = Fraction.from_decimal(numerator)117 self._numerator = value._numerator118 self._denominator = value._denominator119 return self120 121 elif isinstance(numerator, basestring):122 # Handle construction from strings.123 m = _RATIONAL_FORMAT.match(numerator)124 if m is None:125 raise ValueError('Invalid literal for Fraction: %r' %126 numerator)127 numerator = int(m.group('num') or '0')128 denom = m.group('denom')129 if denom:130 denominator = int(denom)131 else:132 denominator = 1133 decimal = m.group('decimal')134 if decimal:135 scale = 10**len(decimal)136 numerator = numerator * scale + int(decimal)137 denominator *= scale138 exp = m.group('exp')139 if exp:140 exp = int(exp)141 if exp >= 0:142 numerator *= 10**exp143 else:144 denominator *= 10**-exp145 if m.group('sign') == '-':146 numerator = -numerator147 148 else:149 raise TypeError(""argument should be a string ""150 ""or a Rational instance"")151 152 elif (isinstance(numerator, Rational) and153 isinstance(denominator, Rational)):154 numerator, denominator = (155 numerator.numerator * denominator.denominator,156 denominator.numerator * numerator.denominator157 )158 else:159 raise TypeError(""both arguments should be ""160 ""Rational instances"")161 162 if denominator == 0:163 raise ZeroDivisionError('Fraction(%s, 0)' % numerator)164 g = gcd(numerator, denominator)165 self._numerator = numerator // g166 self._denominator = denominator // g167 return self168 169 @classmethod170 def from_float(cls, f):171 """"""Converts a finite float to a rational number, exactly.172 173 Beware that Fraction.from_float(0.3) != Fraction(3, 10).174 175 """"""176 if isinstance(f, numbers.Integral):177 return cls(f)178 elif not isinstance(f, float):179 raise TypeError(""%s.from_float() only takes floats, not %r (%s)"" %180 (cls.__name__, f, type(f).__name__))181 if math.isnan(f) or math.isinf(f):182 raise TypeError(""Cannot convert %r to %s."" % (f, cls.__name__))183 return cls(*f.as_integer_ratio())184 185 @classmethod186 def from_decimal(cls, dec):187 """"""Converts a finite Decimal instance to a rational number, exactly.""""""188 from decimal import Decimal189 if isinstance(dec, numbers.Integral):190 dec = Decimal(int(dec))191 elif not isinstance(dec, Decimal):192 raise TypeError(193 ""%s.from_decimal() only takes Decimals, not %r (%s)"" %194 (cls.__name__, dec, type(dec).__name__))195 if not dec.is_finite():196 # Catches infinities and nans.197 raise TypeError(""Cannot convert %s to %s."" % (dec, cls.__name__))198 sign, digits, exp = dec.as_tuple()199 digits = int(''.join(map(str, digits)))200 if sign:201 digits = -digits202 if exp >= 0:203 return cls(digits * 10 ** exp)204 else:205 return cls(digits, 10 ** -exp)206 207 def limit_denominator(self, max_denominator=1000000):208 """"""Closest Fraction to self with denominator at most max_denominator.209 210 >>> Fraction('3.141592653589793').limit_denominator(10)211 Fraction(22, 7)212 >>> Fraction('3.141592653589793').limit_denominator(100)213 Fraction(311, 99)214 >>> Fraction(4321, 8765).limit_denominator(10000)215 Fraction(4321, 8765)216 217 """"""218 # Algorithm notes: For any real number x, define a *best upper219 # approximation* to x to be a rational number p/q such that:220 #221 # (1) p/q >= x, and222 # (2) if p/q > r/s >= x then s > q, for any rational r/s.223 #224 # Define *best lower approximation* similarly. Then it can be225 # proved that a rational number is a best upper or lower226 # approximation to x if, and only if, it is a convergent or227 # semiconvergent of the (unique shortest) continued fraction228 # associated to x.229 #230 # To find a best rational approximation with denominator <= M,231 # we find the best upper and lower approximations with232 # denominator <= M and take whichever of these is closer to x.233 # In the event of a tie, the bound with smaller denominator is234 # chosen. If both denominators are equal (which can happen235 # only when max_denominator == 1 and self is midway between236 # two integers) the lower bound---i.e., the floor of self, is237 # taken.238 239 if max_denominator < 1:240 raise ValueError(""max_denominator should be at least 1"")241 if self._denominator <= max_denominator:242 return Fraction(self)243 244 p0, q0, p1, q1 = 0, 1, 1, 0245 n, d = self._numerator, self._denominator246 while True:247 a = n//d248 q2 = q0+a*q1249 if q2 > max_denominator:250 break251 p0, q0, p1, q1 = p1, q1, p0+a*p1, q2252 n, d = d, n-a*d253 254 k = (max_denominator-q0)//q1255 bound1 = Fraction(p0+k*p1, q0+k*q1)256 bound2 = Fraction(p1, q1)257 if abs(bound2 - self) <= abs(bound1-self):258 return bound2259 else:260 return bound1261 262 @property263 def numerator(a):264 return a._numerator265 266 @property267 def denominator(a):268 return a._denominator269 270 def __repr__(self):271 """"""repr(self)""""""272 return ('Fraction(%s, %s)' % (self._numerator, self._denominator))273 274 def __str__(self):275 """"""str(self)""""""276 if self._denominator == 1:277 return str(self._numerator)278 else:279 return '%s/%s' % (self._numerator, self._denominator)280 281 def _operator_fallbacks(monomorphic_operator, fallback_operator):282 """"""Generates forward and reverse operators given a purely-rational283 operator and a function from the operator module.284 285 Use this like:286 __op__, __rop__ = _operator_fallbacks(just_rational_op, operator.op)287 288 In general, we want to implement the arithmetic operations so289 that mixed-mode operations either call an implementation whose290 author knew about the types of both arguments, or convert both291 to the nearest built in type and do the operation there. In292 Fraction, that means that we define __add__ and __radd__ as:293 294 def __add__(self, other):295 # Both types have numerators/denominator attributes,296 # so do the operation directly297 if isinstance(other, (int, long, Fraction)):298 return Fraction(self.numerator * other.denominator +299 other.numerator * self.denominator,300 self.denominator * other.denominator)301 # float and complex don't have those operations, but we302 # know about those types, so special case them.303 elif isinstance(other, float):304 return float(self) + other305 elif isinstance(other, complex):306 return complex(self) + other307 # Let the other type take over.308 return NotImplemented309 310 def __radd__(self, other):311 # radd handles more types than add because there's312 # nothing left to fall back to.313 if isinstance(other, Rational):314 return Fraction(self.numerator * other.denominator +315 other.numerator * self.denominator,316 self.denominator * other.denominator)317 elif isinstance(other, Real):318 return float(other) + float(self)319 elif isinstance(other, Complex):320 return complex(other) + complex(self)321 return NotImplemented322 323 324 There are 5 different cases for a mixed-type addition on325 Fraction. I'll refer to all of the above code that doesn't326 refer to Fraction, float, or complex as ""boilerplate"". 'r'327 will be an instance of Fraction, which is a subtype of328 Rational (r : Fraction <: Rational), and b : B <:329 Complex. The first three involve 'r + b':330 331 1. If B <: Fraction, int, float, or complex, we handle332 that specially, and all is well.333 2. If Fraction falls back to the boilerplate code, and it334 were to return a value from __add__, we'd miss the335 possibility that B defines a more intelligent __radd__,336 so the boilerplate should return NotImplemented from337 __add__. In particular, we don't handle Rational338 here, even though we could get an exact answer, in case339 the other type wants to do something special.340 3. If B <: Fraction, Python tries B.__radd__ before341 Fraction.__add__. This is ok, because it was342 implemented with knowledge of Fraction, so it can343 handle those instances before delegating to Real or344 Complex.345 346 The next two situations describe 'b + r'. We assume that b347 didn't know about Fraction in its implementation, and that it348 uses similar boilerplate code:349 350 4. If B <: Rational, then __radd_ converts both to the351 builtin rational type (hey look, that's us) and352 proceeds.353 5. Otherwise, __radd__ tries to find the nearest common354 base ABC, and fall back to its builtin type. Since this355 class doesn't subclass a concrete type, there's no356 implementation to fall back to, so we need to try as357 hard as possible to return an actual value, or the user358 will get a TypeError.359 360 """"""361 def forward(a, b):362 if isinstance(b, (int, long, Fraction)):363 return monomorphic_operator(a, b)364 elif isinstance(b, float):365 return fallback_operator(float(a), b)366 elif isinstance(b, complex):367 return fallback_operator(complex(a), b)368 else:369 return NotImplemented370 forward.__name__ = '__' + fallback_operator.__name__ + '__'371 forward.__doc__ = monomorphic_operator.__doc__372 373 def reverse(b, a):374 if isinstance(a, Rational):375 # Includes ints.376 return monomorphic_operator(a, b)377 elif isinstance(a, numbers.Real):378 return fallback_operator(float(a), float(b))379 elif isinstance(a, numbers.Complex):380 return fallback_operator(complex(a), complex(b))381 else:382 return NotImplemented383 reverse.__name__ = '__r' + fallback_operator.__name__ + '__'384 reverse.__doc__ = monomorphic_operator.__doc__385 386 return forward, reverse387 388 def _add(a, b):389 """"""a + b""""""390 return Fraction(a.numerator * b.denominator +391 b.numerator * a.denominator,392 a.denominator * b.denominator)393 394 __add__, __radd__ = _operator_fallbacks(_add, operator.add)395 396 def _sub(a, b):397 """"""a - b""""""398 return Fraction(a.numerator * b.denominator -399 b.numerator * a.denominator,400 a.denominator * b.denominator)401 402 __sub__, __rsub__ = _operator_fallbacks(_sub, operator.sub)403 404 def _mul(a, b):405 """"""a * b""""""406 return Fraction(a.numerator * b.numerator, a.denominator * b.denominator)407 408 __mul__, __rmul__ = _operator_fallbacks(_mul, operator.mul)409 410 def _div(a, b):411 """"""a / b""""""412 return Fraction(a.numerator * b.denominator,413 a.denominator * b.numerator)414 415 __truediv__, __rtruediv__ = _operator_fallbacks(_div, operator.truediv)416 __div__, __rdiv__ = _operator_fallbacks(_div, operator.div)417 418 def __floordiv__(a, b):419 """"""a // b""""""420 # Will be math.floor(a / b) in 3.0.421 div = a / b422 if isinstance(div, Rational):423 # trunc(math.floor(div)) doesn't work if the rational is424 # more precise than a float because the intermediate425 # rounding may cross an integer boundary.426 return div.numerator // div.denominator427 else:428 return math.floor(div)429 430 def __rfloordiv__(b, a):431 """"""a // b""""""432 # Will be math.floor(a / b) in 3.0.433 div = a / b434 if isinstance(div, Rational):435 # trunc(math.floor(div)) doesn't work if the rational is436 # more precise than a float because the intermediate437 # rounding may cross an integer boundary.438 return div.numerator // div.denominator439 else:440 return math.floor(div)441 442 def __mod__(a, b):443 """"""a % b""""""444 div = a // b445 return a - b * div446 447 def __rmod__(b, a):448 """"""a % b""""""449 div = a // b450 return a - b * div451 452 def __pow__(a, b):453 """"""a ** b454 455 If b is not an integer, the result will be a float or complex456 since roots are generally irrational. If b is an integer, the457 result will be rational.458 459 """"""460 if isinstance(b, Rational):461 if b.denominator == 1:462 power = b.numerator463 if power >= 0:464 return Fraction(a._numerator ** power,465 a._denominator ** power)466 else:467 return Fraction(a._denominator ** -power,468 a._numerator ** -power)469 else:470 # A fractional power will generally produce an471 # irrational number.472 return float(a) ** float(b)473 else:474 return float(a) ** b475 476 def __rpow__(b, a):477 """"""a ** b""""""478 if b._denominator == 1 and b._numerator >= 0:479 # If a is an int, keep it that way if possible.480 return a ** b._numerator481 482 if isinstance(a, Rational):483 return Fraction(a.numerator, a.denominator) ** b484 485 if b._denominator == 1:486 return a ** b._numerator487 488 return a ** float(b)489 490 def __pos__(a):491 """"""+a: Coerces a subclass instance to Fraction""""""492 return Fraction(a._numerator, a._denominator)493 494 def __neg__(a):495 """"""-a""""""496 return Fraction(-a._numerator, a._denominator)497 498 def __abs__(a):499 """"""abs(a)""""""500 return Fraction(abs(a._numerator), a._denominator)501 502 def __trunc__(a):503 """"""trunc(a)""""""504 if a._numerator < 0:505 return -(-a._numerator // a._denominator)506 else:507 return a._numerator // a._denominator508 509 def __hash__(self):510 """"""hash(self)511 512 Tricky because values that are exactly representable as a513 float must have the same hash as that float.514 515 """"""516 # XXX since this method is expensive, consider caching the result517 if self._denominator == 1:518 # Get integers right.519 return hash(self._numerator)520 # Expensive check, but definitely correct.521 if self == float(self):522 return hash(float(self))523 else:524 # Use tuple's hash to avoid a high collision rate on525 # simple fractions.526 return hash((self._numerator, self._denominator))527 528 def __eq__(a, b):529 """"""a == b""""""530 if isinstance(b, Rational):531 return (a._numerator == b.numerator and532 a._denominator == b.denominator)533 if isinstance(b, numbers.Complex) and b.imag == 0:534 b = b.real535 if isinstance(b, float):536 if math.isnan(b) or math.isinf(b):537 # comparisons with an infinity or nan should behave in538 # the same way for any finite a, so treat a as zero.539 return 0.0 == b540 else:541 return a == a.from_float(b)542 else:543 # Since a doesn't know how to compare with b, let's give b544 # a chance to compare itself with a.545 return NotImplemented546 547 def _richcmp(self, other, op):548 """"""Helper for comparison operators, for internal use only.549 550 Implement comparison between a Rational instance `self`, and551 either another Rational instance or a float `other`. If552 `other` is not a Rational instance or a float, return553 NotImplemented. `op` should be one of the six standard554 comparison operators.555 556 """"""557 # convert other to a Rational instance where reasonable.558 if isinstance(other, Rational):559 return op(self._numerator * other.denominator,560 self._denominator * other.numerator)561 # comparisons with complex should raise a TypeError, for consistency562 # with int<->complex, float<->complex, and complex<->complex comparisons.563 if isinstance(other, complex):564 raise TypeError(""no ordering relation is defined for complex numbers"")565 if isinstance(other, float):566 if math.isnan(other) or math.isinf(other):567 return op(0.0, other)568 else:569 return op(self, self.from_float(other))570 else:571 return NotImplemented572 573 def __lt__(a, b):574 """"""a < b""""""575 return a._richcmp(b, operator.lt)576 577 def __gt__(a, b):578 """"""a > b""""""579 return a._richcmp(b, operator.gt)580 581 def __le__(a, b):582 """"""a <= b""""""583 return a._richcmp(b, operator.le)584 585 def __ge__(a, b):586 """"""a >= b""""""587 return a._richcmp(b, operator.ge)588 589 def __nonzero__(a):590 """"""a != 0""""""591 return a._numerator != 0592 593 # support for pickling, copy, and deepcopy594 595 def __reduce__(self):596 return (self.__class__, (str(self),))597 598 def __copy__(self):599 if type(self) == Fraction:600 return self # I'm immutable; therefore I am my own clone601 return self.__class__(self._numerator, self._denominator)602 603 def __deepcopy__(self, memo):604 if type(self) == Fraction:605 return self # My components are also immutable606 return self.__class__(self._numerator, self._denominator)607 608","Python"
609"Multi-scale modeling","vishalsubbiah/Grain-Boundary-Energies-LAMMPS","Code and Scripts/Python and Lammps/FullStackAll/FullStack101010/Experiments/Fe/300/WriteLammps.py",".py","5580","138","import shutil610 611def write_lammps(M1,M2,potential,element,temp,tempDamp,theta,lattice,minEng,typeAtom):612 613 x1=[int(M1[0][0]),int(M1[0][1]),int(M1[0][2])]614 y1=[int(M1[1][0]),int(M1[1][1]),int(M1[1][2])]615 z1=[int(M1[2][0]),int(M1[2][1]),int(M1[2][2])]616 617 x2=[int(M2[0][0]),int(M2[0][1]),int(M2[0][2])]618 y2=[int(M2[1][0]),int(M2[1][1]),int(M2[1][2])]619 z2=[int(M2[2][0]),int(M2[2][1]),int(M2[2][2])]620 621 shutil.copy('../../Potentials/'+str(potential),'.')622 623 LammpsInFile=open('GB.in','w')624 625 ContentComment=""""""# LAMMPS Input File for Grain Boundaries 626# Mark Tschopp, Dec2009 627# This file will generate a single Sigma5(310) STGB \n""""""628 629 ContentInitSim=""""""\n # ---------- Initialize Simulation --------------------- 630clear 631units metal 632dimension 3 633boundary p p p 634atom_style atomic \n""""""635 636 ContentAtomStruct1=""""""\n # ---------- Create Atomistic Structure --------------------- 637lattice """"""+str(typeAtom)+"""""" """"""+str(lattice)+"""""" 638region whole block -10 10 -10 10 -10 10 units lattice 639create_box 2 whole 640region upper block INF INF 0 10 INF INF units lattice 641 642lattice """"""+str(typeAtom)+"""""" """"""+str(lattice)+"""""" orient x """""" +str(x1[0])+"" "" + str(x1[1]) + "" "" + str(x1[2]) +"""""" orient y """"""+ str(y1[0])+"" "" + str(y1[1]) + "" "" + str(y1[2]) +"""""" orient z """""" +str(z1[0])+"" "" + str(z1[1]) + "" "" + str(z1[2]) +""""""\n"""""" 643 644 645 ContentAtomStruct2=""""""\n create_atoms 1 region upper 646region lower block INF INF -10 0.000000 INF INF units lattice 647lattice """"""+str(typeAtom)+"""""" """"""+str(lattice)+"""""" orient x """""" +str(x2[0])+"" "" + str(x2[1]) + "" "" + str(x2[2]) +"""""" orient y """"""+ str(y2[0])+"" "" + str(y2[1]) + "" "" + str(y2[2]) +"""""" orient z """""" +str(z2[0])+"" "" + str(z2[1]) + "" "" + str(z2[2]) +""""""\n""""""648 649 ContentAtomStruct3=""""""\n create_atoms 2 region lower 650group upper type 1 651group lower type 2 652replicate 1 1 1 \n""""""653 654 ContentInterPot=""""""\n # ---------- Define Interatomic Potential --------------------- 655pair_style eam/alloy 656pair_coeff * * """"""+str(potential)+"" ""+str(element)+ "" ""+str(element) +""""""657neighbor 2.0 bin 658neigh_modify delay 10 check yes \n""""""659 660 ContentDisplace=""""""\n # ---------- Displace atoms and delete overlapping atoms --------------------- 661displace_atoms upper move 0 0 0 units lattice 662delete_atoms overlap 0.35 lower upper \n""""""663 664 ContentSettings=""""""\n # ---------- Define Settings --------------------- 665compute csym all centro/atom """"""+str(typeAtom)+""""""666compute eng all pe/atom 667compute eatoms all reduce sum c_eng 668 669#write_dump all custom dump.img_orig.cfg mass type xs ys zs \n""""""670 671 ContentEquil=""""""\n # ----------- EQUILIBRATION --------------------672write_dump all custom dump.img_pre_equil_""""""+str(theta)+"""""".cfg mass type xs ys zs673reset_timestep 0674timestep 0.001675velocity all create """"""+str(temp)+"""""" 12345 mom yes rot no676fix 2 all npt temp """"""+ str(temp) + "" "" + str(temp)+ "" "" + str(tempDamp)+"""""" iso 0 0 1 drag 1 677 678thermo 1000 679thermo_style custom step temp pe lx ly lz press pxx pyy pzz c_eatoms 680dump 5 all cfg 1000 dump.equal_*.cfg mass type xs ys zs c_csym c_eng fx fy fz681dump_modify 5 element """"""+str(element)+"" ""+str(element) + """"""682 683# Run for at least 10 picosecond (assuming 1 fs timestep)684#write_dump all custom dump.img_post_equil.cfg mass type xs ys zs685 686run 100000687unfix 2\n""""""688 689 ContentMin1=""""""\n # ---------- Run Minimization --------------------- 690write_dump all custom dump.img_post_equil_""""""+str(theta)+"""""".cfg mass type xs ys zs 691reset_timestep 0 692thermo 10 693thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms 694dump 1 all cfg 25 dump.sig5_minimization_*.cfg mass type xs ys zs c_csym c_eng fx fy fz695dump_modify 1 element """"""+str(element)+"" ""+str(element)+""""""696min_style cg 697minimize 1e-15 1e-15 5000 5000 698undump 1 \n""""""699 700 ContentMin2=""""""\n # ---------- Run Minimization 2--------------------- 701# Now allow the box to expand/contract perpendicular to the grain boundary702reset_timestep 0 703thermo 10 704thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms 705fix 1 all box/relax y 0 vmax 0.001706min_style cg 707minimize 1e-15 1e-15 5000 5000 \n""""""708 709 ContentGBEng=""""""\n # ---------- Calculate GB Energy --------------------- 710variable minimumenergy equal """"""+str(minEng)+""""""711variable esum equal ""v_minimumenergy * count(all)"" 712variable xseng equal ""c_eatoms - (v_minimumenergy * count(all))"" 713variable gbarea equal ""lx * lz * 2"" 714variable gbe equal ""(c_eatoms - (v_minimumenergy * count(all)))/v_gbarea"" 715variable gbemJm2 equal ${gbe}*16021.7733 716variable gbernd equal round(${gbemJm2}) 717print ""GB energy is ${gbemJm2} mJ/m^2"" \n""""""718 719 ContentDumpData=""""""\n # ---------- Dump data into Data file ------------- 720reset_timestep 0 721dump 1 all cfg 1000 dump.al_sig5_310_*.cfg mass type xs ys zs c_csym c_eng fx fy fz722dump_modify 1 element """"""+str(element)+"" ""+str(element)+""""""723minimize 1e-15 1e-15 5000 5000724undump 1725 726write_restart restart.al_sig5_310_stgb727write_dump all custom dump.img_post_minimize_""""""+str(theta)+"""""".cfg mass type xs ys zs728print ""All done"" \n""""""729 730 LammpsInFile.write(ContentComment)731 LammpsInFile.write(ContentInitSim)732 LammpsInFile.write(ContentAtomStruct1)733 LammpsInFile.write(ContentAtomStruct2)734 LammpsInFile.write(ContentAtomStruct3)735 LammpsInFile.write(ContentInterPot)736 LammpsInFile.write(ContentDisplace)737 LammpsInFile.write(ContentSettings)738 if(temp != 0):739 LammpsInFile.write(ContentEquil)740 LammpsInFile.write(ContentMin1)741 LammpsInFile.write(ContentMin2)742 LammpsInFile.write(ContentGBEng)743 LammpsInFile.write(ContentDumpData)744 745 LammpsInFile.close()746","Python"
747"Multi-scale modeling","vishalsubbiah/Grain-Boundary-Energies-LAMMPS","Code and Scripts/Python and Lammps/FullStackAll/FullStack101010/Experiments/Fe/300/Quaternion.py",".py","12278","415","""""""748Quaternion provides a class for manipulating quaternion objects. This class provides:749 750 - a convenient constructor to convert to/from Euler Angles (RA,Dec,Roll) 751 to/from quaternions752 - class methods to multiply and divide quaternions 753""""""754 755__copyright__ = """"""756Copyright (c) 2009, Smithsonian Astrophysical Observatory757All rights reserved.758 759Redistribution and use in source and binary forms, with or without760modification, are permitted provided that the following conditions are met:761 * Redistributions of source code must retain the above copyright762 notice, this list of conditions and the following disclaimer.763 * Redistributions in binary form must reproduce the above copyright764 notice, this list of conditions and the following disclaimer in the765 documentation and/or other materials provided with the distribution.766 * Neither the name of the <organization> nor the767 names of its contributors may be used to endorse or promote products768 derived from this software without specific prior written permission.769 770THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ""AS IS"" AND771ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED772WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE773DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY774DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES775(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;776LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND777ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT778(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS779SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.780""""""781 782 783 784 785import numpy as np786from math import cos, sin, radians, degrees, atan2, sqrt787 788class Quat(object):789 """"""790 Quaternion class791 792 Example usage::793 794 >>> from Quaternion import Quat795 >>> quat = Quat((12,45,45))796 >>> quat.ra, quat.dec, quat.roll797 (12, 45, 45)798 >>> quat.q799 array([ 0.38857298, -0.3146602 , 0.23486498, 0.8335697 ])800 >>> q2 = Quat([ 0.38857298, -0.3146602 , 0.23486498, 0.8335697])801 >>> q2.ra802 11.999999315925008803 804 805 Multiplication and division operators are overloaded for the class to 806 perform appropriate quaternion multiplication and division.807 808 Example usage::809 810 >>> q1 = Quat((20,30,40))811 >>> q2 = Quat((30,40,50))812 >>> q = q1 / q2813 814 Performs the operation as q1 * inverse q2815 816 Example usage::817 818 >>> q1 = Quat((20,30,40))819 >>> q2 = Quat((30,40,50))820 >>> q = q1 * q2821 822 823 :param attitude: initialization attitude for quat824 825 ``attitude`` may be:826 * another Quat827 * a 4 element array (expects x,y,z,w quat form)828 * a 3 element array (expects ra,dec,roll in degrees)829 * a 3x3 transform/rotation matrix830 831 """"""832 def __init__(self, attitude):833 self._q = None834 self._equatorial = None835 self._T = None836 # checks to see if we've been passed a Quat 837 if isinstance(attitude, Quat):838 self._set_q(attitude.q)839 else:840 # make it an array and check to see if it is a supported shape841 attitude = np.array(attitude)842 if len(attitude) == 4:843 self._set_q(attitude)844 elif attitude.shape == (3,3):845 self._set_transform(attitude)846 elif attitude.shape == (3,):847 self._set_equatorial(attitude)848 else:849 raise TypeError(""attitude is not one of possible types (3 or 4 elements, Quat, or 3x3 matrix)"")850 851 852 def _set_q(self, q):853 """"""854 Set the value of the 4 element quaternion vector 855 856 :param q: list or array of normalized quaternion elements857 """"""858 q = np.array(q)859 if abs(np.sum(q**2) - 1.0) > 1e-6:860 raise ValueError('Quaternion must be normalized so sum(q**2) == 1; use Quaternion.normalize')861 self._q = (q if q[3] > 0 else -q)862 # Erase internal values of other representations863 self._equatorial = None864 self._T = None865 866 def _get_q(self):867 """"""868 Retrieve 4-vector of quaternion elements in [x, y, z, w] form869 870 :rtype: numpy array871 872 """"""873 if self._q is None:874 # Figure out q from available values, doing nothing others are not defined875 if self._equatorial is not None:876 self._q = self._equatorial2quat()877 elif self._T is not None:878 self._q = self._transform2quat()879 return self._q880 881 # use property to make this get/set automatic882 q = property(_get_q, _set_q)883 884 def _set_equatorial(self, equatorial):885 """"""Set the value of the 3 element equatorial coordinate list [RA,Dec,Roll]886 expects values in degrees887 bounds are not checked888 889 :param equatorial: list or array [ RA, Dec, Roll] in degrees890 891 """"""892 att = np.array(equatorial)893 ra, dec, roll = att894 self._ra0 = ra895 if ( ra > 180 ):896 self._ra0 = ra - 360897 self._roll0 = roll898 if ( roll > 180):899 self._roll0 = roll - 360900 self._equatorial = att901 902 def _get_equatorial(self):903 """"""Retrieve [RA, Dec, Roll]904 905 :rtype: numpy array906 """"""907 if self._equatorial is None:908 if self._q is not None:909 self._equatorial = self._quat2equatorial()910 elif self._T is not None:911 self._q = self._transform2quat()912 self._equatorial = self._quat2equatorial()913 return self._equatorial914 915 equatorial = property(_get_equatorial,_set_equatorial)916 917 def _get_ra(self):918 """"""Retrieve RA term from equatorial system in degrees""""""919 return self.equatorial[0]920 921 def _get_dec(self):922 """"""Retrieve Dec term from equatorial system in degrees""""""923 return self.equatorial[1]924 925 def _get_roll(self):926 """"""Retrieve Roll term from equatorial system in degrees""""""927 return self.equatorial[2]928 929 ra = property(_get_ra)930 dec = property(_get_dec)931 roll = property(_get_roll)932 933 def _set_transform(self, T):934 """"""935 Set the value of the 3x3 rotation/transform matrix936 937 :param T: 3x3 array/numpy array938 """"""939 transform = np.array(T)940 self._T = transform941 942 def _get_transform(self):943 """"""944 Retrieve the value of the 3x3 rotation/transform matrix945 946 :returns: 3x3 rotation/transform matrix947 :rtype: numpy array948 949 """"""950 if self._T is None:951 if self._q is not None:952 self._T = self._quat2transform()953 elif self._equatorial is not None:954 self._T = self._equatorial2transform()955 return self._T956 957 transform = property(_get_transform, _set_transform)958 959 def _quat2equatorial(self):960 """"""961 Determine Right Ascension, Declination, and Roll for the object quaternion962 963 :returns: RA, Dec, Roll964 :rtype: numpy array [ra,dec,roll]965 """"""966 967 q = self.q968 q2 = self.q**2969 970 ## calculate direction cosine matrix elements from $quaternions971 xa = q2[0] - q2[1] - q2[2] + q2[3] 972 xb = 2 * (q[0] * q[1] + q[2] * q[3]) 973 xn = 2 * (q[0] * q[2] - q[1] * q[3]) 974 yn = 2 * (q[1] * q[2] + q[0] * q[3]) 975 zn = q2[3] + q2[2] - q2[0] - q2[1] 976 977 ##; calculate RA, Dec, Roll from cosine matrix elements978 ra = degrees(atan2(xb , xa)) ;979 dec = degrees(atan2(xn , sqrt(1 - xn**2)));980 roll = degrees(atan2(yn , zn)) ;981 if ( ra < 0 ):982 ra += 360983 if ( roll < 0 ):984 roll += 360985 986 return np.array([ra, dec, roll])987 988 989 def _quat2transform(self):990 """"""991 Transform a unit quaternion into its corresponding rotation matrix (to992 be applied on the right side).993 994 :returns: transform matrix995 :rtype: numpy array996 997 """"""998 x, y, z, w = self.q999 xx2 = 2 * x * x1000 yy2 = 2 * y * y1001 zz2 = 2 * z * z1002 xy2 = 2 * x * y1003 wz2 = 2 * w * z1004 zx2 = 2 * z * x1005 wy2 = 2 * w * y1006 yz2 = 2 * y * z1007 wx2 = 2 * w * x1008 1009 rmat = np.empty((3, 3), float)1010 rmat[0,0] = 1. - yy2 - zz21011 rmat[0,1] = xy2 - wz21012 rmat[0,2] = zx2 + wy21013 rmat[1,0] = xy2 + wz21014 rmat[1,1] = 1. - xx2 - zz21015 rmat[1,2] = yz2 - wx21016 rmat[2,0] = zx2 - wy21017 rmat[2,1] = yz2 + wx21018 rmat[2,2] = 1. - xx2 - yy21019 1020 return rmat1021 1022 def _equatorial2quat( self ):1023 """"""Dummy method to return return quat. 1024 1025 :returns: quaternion1026 :rtype: Quat1027 1028 """"""1029 return self._transform2quat()1030 1031 def _equatorial2transform( self ):1032 """"""Construct the transform/rotation matrix from RA,Dec,Roll1033 1034 :returns: transform matrix1035 :rtype: 3x3 numpy array1036 1037 """"""1038 ra = radians(self._get_ra())1039 dec = radians(self._get_dec())1040 roll = radians(self._get_roll())1041 ca = cos(ra)1042 sa = sin(ra)1043 cd = cos(dec)1044 sd = sin(dec)1045 cr = cos(roll)1046 sr = sin(roll)1047 1048 # This is the transpose of the transformation matrix (related to translation1049 # of original perl code1050 rmat = np.array([[ca * cd, sa * cd, sd ],1051 [-ca * sd * sr - sa * cr, -sa * sd * sr + ca * cr, cd * sr],1052 [-ca * sd * cr + sa * sr, -sa * sd * cr - ca * sr, cd * cr]])1053 1054 return rmat.transpose()1055 1056 def _transform2quat( self ):1057 """"""Construct quaternion from the transform/rotation matrix 1058 1059 :returns: quaternion formed from transform matrix1060 :rtype: numpy array1061 """"""1062 1063 # Code was copied from perl PDL code that uses backwards index ordering1064 T = self.transform.transpose() 1065 den = np.array([ 1.0 + T[0,0] - T[1,1] - T[2,2],1066 1.0 - T[0,0] + T[1,1] - T[2,2],1067 1.0 - T[0,0] - T[1,1] + T[2,2],1068 1.0 + T[0,0] + T[1,1] + T[2,2]])1069 1070 max_idx = np.flatnonzero(den == max(den))[0]1071 1072 q = np.zeros(4)1073 q[max_idx] = 0.5 * sqrt(max(den))1074 denom = 4.0 * q[max_idx]1075 if (max_idx == 0):1076 q[1] = (T[1,0] + T[0,1]) / denom 1077 q[2] = (T[2,0] + T[0,2]) / denom 1078 q[3] = -(T[2,1] - T[1,2]) / denom 1079 if (max_idx == 1):1080 q[0] = (T[1,0] + T[0,1]) / denom 1081 q[2] = (T[2,1] + T[1,2]) / denom 1082 q[3] = -(T[0,2] - T[2,0]) / denom 1083 if (max_idx == 2):1084 q[0] = (T[2,0] + T[0,2]) / denom 1085 q[1] = (T[2,1] + T[1,2]) / denom 1086 q[3] = -(T[1,0] - T[0,1]) / denom 1087 if (max_idx == 3):1088 q[0] = -(T[2,1] - T[1,2]) / denom 1089 q[1] = -(T[0,2] - T[2,0]) / denom 1090 q[2] = -(T[1,0] - T[0,1]) / denom 1091 1092 return q1093 1094 1095 def __div__(self, quat2):1096 """"""1097 Divide one quaternion by another.1098 1099 Example usage::1100 1101 >>> q1 = Quat((20,30,40))1102 >>> q2 = Quat((30,40,50))1103 >>> q = q1 / q21104 1105 Performs the operation as q1 * inverse q21106 1107 :returns: product q1 * inverse q21108 :rtype: Quat1109 1110 """"""1111 return self * quat2.inv()1112 1113 1114 def __mul__(self, quat2):1115 """"""1116 Multiply quaternion by another.1117 1118 Example usage::1119 1120 >>> q1 = Quat((20,30,40))1121 >>> q2 = Quat((30,40,50))1122 >>> (q1 * q2).equatorial1123 array([ 349.73395729, 76.25393056, 127.61636727])1124 1125 :returns: product q1 * q21126 :rtype: Quat1127 1128 """"""1129 q1 = self.q1130 q2 = quat2.q1131 mult = np.zeros(4)1132 mult[0] = q1[3]*q2[0] - q1[2]*q2[1] + q1[1]*q2[2] + q1[0]*q2[3]1133 mult[1] = q1[2]*q2[0] + q1[3]*q2[1] - q1[0]*q2[2] + q1[1]*q2[3]1134 mult[2] = -q1[1]*q2[0] + q1[0]*q2[1] + q1[3]*q2[2] + q1[2]*q2[3]1135 mult[3] = -q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] + q1[3]*q2[3]1136 return Quat(mult)1137 1138 def inv(self):1139 """"""1140 Invert the quaternion 1141 1142 :returns: inverted quaternion1143 :rtype: Quat1144 """"""1145 return Quat([self.q[0], self.q[1], self.q[2], -self.q[3]])1146 1147 1148def normalize(array):1149 """""" 1150 Normalize a 4 element array/list/numpy.array for use as a quaternion1151 1152 :param quat_array: 4 element list/array1153 :returns: normalized array1154 :rtype: numpy array1155 1156 """"""1157 quat = np.array(array)1158 return quat / np.sqrt(np.dot(quat, quat))1159 1160 1161","Python"
1162"Multi-scale modeling","vishalsubbiah/Grain-Boundary-Energies-LAMMPS","Code and Scripts/Python and Lammps/FullStackAll/FullStack101010/Experiments/Fe/300/ReadDump.py",".py","617","20","nData=[]1163AData=[]1164for i in range(0,1):1165 filename=""dump.img_post_minimize_""+str(i)+"".cfg""1166 with open(filename, 'r') as fileData:1167 for j in range(0,7):1168 if(j==3): 1169 nData.append(int(fileData.readline().strip()))1170 elif(j==6):1171 A=fileData.readline().strip().split(' ')1172 AData.append(float(A[0])*float(A[1])*-1)1173 else:1174 fileData.readline()1175 1176filename = ""AreaAtoms.txt""1177 1178with open(filename, 'w') as fileData:1179 for i in range(0,1):1180 fileData.writelines(str(nData[i]) + "" "" + str(AData[i]) + ""\n"")1181","Python"
1182"Multi-scale modeling","vishalsubbiah/Grain-Boundary-Energies-LAMMPS","Code and Scripts/Python and Lammps/FullStackAll/FullStack101010/Experiments/Fe/300/ReadLog.py",".py","251","13","from re import findall1183 1184def read_log():1185 log=open('log.lammps','r')1186 for i in log.readlines():1187 t = i.find('GB energy')1188 if(t == 0):1189 Energy=i.split(' ')1190 GBEnergy = float(Energy[3])1191 return GBEnergy1192 1193 1194","Python"
1195"Multi-scale modeling","vishalsubbiah/Grain-Boundary-Energies-LAMMPS","Code and Scripts/Python and Lammps/FullStackAll/FullStack101010/Experiments/Fe/300/RunLammps.py",".py","1100","44","from Quaternion import Quat, normalize1196from WriteLammps import write_lammps1197from QuantFunc import QuatMat, find_orient, scale1198from ReadLog import read_log1199#from GenerateAtoms import GenerateAtoms1200import os