Ask Your Question

RyanG's profile - activity

2020-09-19 17:24:51 -0600 asked a question Unable to solve basic equation, what am I doing wrong?

I'm trying to solve a simple exponential equation, but Sage freezes on every attempt until I manually stop the process:

 sage: var('E_j','T'); k = 0.695201368776488;
 sage: Boltz(E_j,T) = e^(-(1/(k*T))*E_j)
 sage: f = Boltz(300,T)==0.15
 sage: f
 e^(-431.529645184649/T) == 0.150000000000000

neither

f.solve(T)

or

f.solve(T,to_poly_solve=True)

Provide a solution and force me to interrupt the console with ctr-C

2020-06-28 10:53:14 -0600 received badge  Popular Question (source)
2020-06-06 07:21:58 -0600 asked a question Creating a supercell from a unit cell

I'm trying to figure out how to create an arbitrary three dimensional supercell from a 6 coordinate unit cell using Sage. I'm hoping there is someone who may be able to suggest a few functions that could help me solve this problem. Currently I'm approaching it like this but am having trouble making it scalable to other dimensions.

x = 0.3
h=1/2
q = 1/4

#Rutile has six atoms in its primitive cell
#Titanium atoms (0.0, 0.0, 0.0), (h, h, h)
#Oxygen atoms ((x, x, 0.0), (1-x, 1-x, 0.0), (-x+h, x+h, h), (x+h, -x+h, h)
#these coordinates are for a cell that is doubled in x and y'

#Coordinates for a 2 x 2 x 1 supercell

#layer 1 Ti: 

(0.0, 0.0, 0.0), (h, 0.0, 0.0), (0.0, h, 0.0), ('Ti', h, h, 0.0)

#layer 1 O: 

(xh, xh , 0.0),  (xh, xh+h, 0.0),  (xh+h, xh, 0.0),  (xh+h, xh+h, 0.0), (1-xh, 1-xh, 0.0)

(1-xh, 1-xh-h, 0.0), (1-xh-h, 1-xh, 0.0), (1-xh-h, 1-xh-h, 0.0)

#layer 2 Ti

(q, q, h), (q, q+h, h), (q+h, q, h), (q+h, q+h, h)

#layer 2 O

((-x+h)/2.0, (x+h)/2.0, h), ((-x+h)/2.0, h+(x+h)/2.0, h), (h+(-x+h)/2.0, (x+h)/2.0, h), (h+(-x+h)/2.0,h+(x+h)/2.0, h)

((x+h)/2.0,(-x+h)/2.0, h),  (x+h)/2.0, h+(-x+h)/2.0, h), (h+(x+h)/2.0, (-x+h)/2.0, h), ( h+(x+h)/2.0, h+(-x+h)/2.0, h)
2020-05-22 17:20:30 -0600 received badge  Supporter (source)
2020-05-22 14:28:58 -0600 commented answer plotting regions in 3D

This was very helpful for me. Thank you.

2020-04-25 19:15:23 -0600 received badge  Commentator
2020-04-25 19:15:23 -0600 commented question How to save 3d plot as rotatable html file

I included another tip as an answer. I was having trouble pasting the code as a comment

2020-04-25 19:14:44 -0600 answered a question How to save 3d plot as rotatable html file

Running the code with a Sage notebook appears to have fixed the issue. Note that if you have want to view the HTML files on a different computer or host them on a website etc, make sure you update the the following code HTML in your three.js file:

Original file HTML (lines 21-24)

<body>
<script src="/Applications/SageMath/local/share/threejs/build/three.min.js"></script>
<script src="/Applications/SageMath/local/share/threejs/examples/js/controls/OrbitControls.js"></script>
<script>

Updated HTML:

<body>
<script src="https://cdn.jsdelivr.net/gh/mrdoob/three.js@r110/build/three.min.js"></script>
<script src="https://cdn.jsdelivr.net/gh/mrdoob/three.js@r110/examples/js/controls/OrbitControls.js"></script>
<script>
2020-04-25 19:03:47 -0600 commented question How to save 3d plot as rotatable html file

I actually coded my own framework for viewing/manipulating molecular structures with Python/Sage. I've never used Jmol at all but I wouldn't be opposed to it. I'm glad it worked for you!

2020-04-24 22:37:32 -0600 commented question How to save 3d plot as rotatable html file

I'm running Sage 9.0 on my MacBook and much of my work involves producing 3D crystal structures (three.js) files which can be saved as HTML (an example being the above link). Maybe update Sage to 9.0 with Homebrew on your Mac?

2020-04-20 21:13:45 -0600 commented question How to save 3d plot as rotatable html file

Could you paste the code associated with the 3D image directly into a Sage terminal? That would launch a graphics window like this:

http://chem302.xyz/wp-content/uploads...

The icon in the bottom right hand corner allows you to save an image directly as html.

2020-04-20 21:12:37 -0600 answered a question How to save 3d plot as rotatable html file

Could you paste the code associated with the 3D image directly into a Sage terminal? That would launch a graphics window like this:

http://chem302.xyz/wp-content/uploads...

The icon in the bottom right hand corner allows you to save an image directly as html.

2020-04-17 21:51:31 -0600 commented question Using implicit_plot3d to create a plane between three points

Thank you, polygon3d worked perfectly. HyperplaneArrangements looks very useful as well.

2020-04-16 13:18:03 -0600 asked a question Using implicit_plot3d to create a plane between three points

Hello all. I'm an undergrad at the beginning phases of an applied crystallography/group theory research project and have been using Sage to make some simple visualization/symmetry analysis scripts. What I'm currently trying to do is create a plane for each face of certain high symmetry groups. Here's an example:

http://chem302.xyz/wp-content/uploads...

How would I create a plane to represent each face of the tetrahedron? In the code below S[1] through S[4] represent the different coordinates in space that make up the substituent molecules/vertices.

def tetra(XY4,q):

    A = { 0: (0,0,0) }

    S = {3 : (2**(1/2), 1, 0), 2 : (-2**(1/2), 1, 0), 1: (0, -1, 2**(1/2)), 4: (0, -1, -2**(1/2))}

    C2R = { 1 : (0,2.5,0), 2: (0,-2.5,0) }

    C3R = { 1 : (0, -2*.75, 2*2**(1/2)*.75), 2 : (0, 2*.75, -2*2**(1/2)*.75) }

    C2R_pair = { 1 : [ C2R[1],C2R[2] ] }

    C3R_pair = { 2 : [ C3R[1],C3R[2] ] }

    Axis_label = { (0, (-2*.75)-0.10, (2*2**(1/2)*.75)+0.10): 'C3', (0,2.75,0): 'C2' }

    Substituents = { S[1]: 'H', S[2]: 'H', S[3]: 'H', S[4]: 'H' }

    Additions = { S[1]: 'H', S[2]: 'H', S[3]: 'H', S[4]: 'H' }

    Central_Atom = { A[0] : 'N' }

    Sub_Bonds = { 'A0-S1' : [A[0],S[1]], 'A0-S2' : [A[0],S[2]], 'A0-S3' : [A[0],S[3]], 'A0-S4': [A[0],S[4]] }

    Sym_outline = {'S2-S3': [S[2],S[3]], 'S3-S4' : [S[3],S[4]], 'S2-S4': [S[2],S[4]], 'S1-S2': [S[1],S[2]], 'S1-S3': [S[1],S[3]], 'S1-S4': [S[1],S[4]]}

    Sym_Scale = { 'yz' : [(-1.5,1.5),(-2**(1/2),2**(1/2)),(-2**(1/2),2**(1/2))], 'xy': [(-1.5,1.5),(-2,2),(-2,2)] }

So if I want to create a plane between:

S[1] = (0, -1, 2**(1/2), S[2] = (-2**(1/2), 1, 0), S[3] = (2**(1/2), 1, 0)

Would I use implicit_plot3d? I tried to use list_plot3d but it did not cover the entire face.

2020-03-18 03:43:46 -0600 received badge  Self-Learner (source)
2020-03-18 03:43:46 -0600 received badge  Teacher (source)
2020-03-16 22:19:36 -0600 commented question Is it possible to determine if a line segment is present between two points on a 3d graphic object?

I have fixed the issue and updated the code. Thank you again.

2020-03-16 17:31:35 -0600 commented question Is it possible to determine if a line segment is present between two points on a 3d graphic object?

I think this will be a much better approach and I imagine any answer to that question would have taken me down a much more complicated rabbit hole.

2020-03-16 17:10:00 -0600 answered a question Is it possible to determine if a line segment is present between two points on a 3d graphic object?

Using John Palmieri's suggestion I created a default dictionary of single bonds and an empty dictionary for pi bonds:

When a pi bond pair is added to Pi_Dict, the pairing from Sigma_Dict is deleted

I was not aware that lists could be entered into dictionaries and this will make my work much easier.

Updated code:

import sys
from sage.plot.plot3d.shapes import LineSegment, Sphere
from sage.plot.colors import rainbow
from sage.plot.plot3d.implicit_plot3d import implicit_plot3d
from sage.plot.plot3d.shapes2 import text3d
from sage.repl.rich_output.pretty_print import show

Radii = {'O': 0.73, 'N': 0.75, 'C': 0.77, 'He': 0.32, 'H': 0.37, 'S': 1.02, 'Cl': 0.99, 
    'F': 0.71, 'Xe': 1.30, 'Si': 1.11, 'B': 0.82, 'P': 1.06, 'Br': 1.14}

Valency = {'O': 6, 'N': 5, 'C': 4, 'He': 0, 'H': 1, 'S': 6, 'Cl': 7, 
    'F': 7, 'Xe': 8, 'Si': 4, 'B': 3, 'P': 5, 'Br': 7}

Colors = {'O': 'red', 'N': 'green', 'C': 'black', 'He': 'cyan', 'H': 'white', 'S': 'yellow',
    'Si': 'purple', 'Xe': 'pink', 'F': 'blue', 'Cl': 'green', 'B': 'pink', 'P': 'orange',
    'Br': 'red'}

X = int(input('Enter x axis shift or press enter for default: '))

def distance(A,B):
    D = float(((B[0]-A[0])**2 + (B[1]-A[1])**2 + (B[2]-A[2])**2)**(1/2))
    return D
def Normalize(x):
    x = (float(x)/0.77)*0.20
    return x

O = (X,0,0)

A2 = (X, 1, 0.559)
A3 = (X, 1, -0.559)
A5 = (X, -1, -0.559)
A6 = (X, -1, 0.559)
A1 = (X, 0, distance(A2,A3))
A4 = (X, 0,-1*distance(A2,A3))

Sigma_Dict = { 'A1-A2' : [A1,A2], 'A2-A3' : [A2,A3], 'A3-A4' : [A3,A4], 'A4-A5' : [A4,A5], 'A5-A6' : [A5,A6], 'A1-A6' : [A1, A6] }

Pi_Dict = { 'A1-A2' : [A1,A2], 'A5-A6' : [A5,A6], 'A3-A4': [A3,A4] }

Atomic_Position = { A1: 'Br', A2: 'C', A3: 'Cl', A4: 'C', A5: 'B', A6: 'C' }

def bond_split(AxAy):
    x = AxAy[0]
    y = AxAy[1]
    return(x,y)

def pi_bond(Cx,Cy):
    X = 0.05
    if Cx[1] == Cy[1]:
        Ca = Cx[1]-X
        Cb = Cx[1]+X
        Cc = Cy[1]-X
        Cd = Cy[1]+X
        Cx1 = (Cx[0], Ca, Cx[2])
        Cx2 = (Cx[0], Cb, Cx[2])
        Cy1 = (Cy[0], Cc, Cy[2])
        Cy2 = (Cy[0], Cd, Cy[2])
        D = LineSegment(Cx1, Cy1, 1, color='white', axes=False, frame=False)
        D += LineSegment(Cx2, Cy2, 1, color='white', axes=False, frame=False)
        return (D)
    else:
        Ca = Cx[2]-X
        Cb = Cx[2]+X
        Cc = Cy[2]-X
        Cd = Cy[2]+X
        Cx1 = (Cx[0], Cx[1], Ca)
        Cx2 = (Cx[0], Cx[1], Cb)
        Cy1 = (Cy[0], Cy[1], Cc)
        Cy2 = (Cy[0], Cy[1], Cd)
        D = LineSegment(Cx1, Cy1, 1, color='white', axes=False, frame=False)
        D += LineSegment(Cx2, Cy2, 1, color='white', axes=False, frame=False)
        return (D)

Benzene = Sphere(.000001, color='white').translate(O)

for atom in Atomic_Position:
    Benzene += Sphere(Normalize(Radii[Atomic_Position[atom]]), color=Colors[Atomic_Position[atom]]).translate(atom)

for pi in Pi_Dict:
    if pi in Sigma_Dict:
        del Sigma_Dict[pi]

for pi in Pi_Dict:
    Benzene += pi_bond(bond_split(Pi_Dict[pi])[0],bond_split(Pi_Dict[pi])[1])

for bond in Sigma_Dict:
    Benzene += LineSegment(bond_split(Sigma_Dict[bond])[0],bond_split(Sigma_Dict[bond])[1], 1, color= 'white' )

show(Benzene)
2020-03-16 15:59:10 -0600 commented question Is it possible to determine if a line segment is present between two points on a 3d graphic object?

Am I understanding correctly that you're suggesting I create a default position dictionary that would like something like:

graph = { A1 : [A2, A6], A2 : [A1, A3], A3 : [A2, A4], A4: [A3, A5], A5: [A4, A6], A6: [A5, A1] }

This way when a user adds a pi bond pair ie A1-A2 I can delete the pairing from the default dictionary and add them to a second dictionary for double bonds or something to that effect. Much more elegant. Thank you for the suggestion.

2020-03-16 14:22:28 -0600 asked a question Is it possible to determine if a line segment is present between two points on a 3d graphic object?

Hello everyone. I'm working on a program which plots 3D molecular structures and associated point group symmetries. I've had some issues figuring out how to fill in the remaining single bonds of a six member ring (aromatic class) after the user dictates where double bonds should be placed. These aromatics are basically 6 member cycle graphs for the purposes of my problem. Is there any way to evaluate two coordinates (atoms) on a 3d graphic object (molecule) and check if a line segment (bond) has already been added between them? By necessity any pi (double) bond should only be connected to the adjacent, (non pi-bonded) atom. So for example if there was a pi bond present between Atoms 1 and 2 of a ring then we should have single bonds between A6-A1 and A2-A3. The code below is how I am currently approaching that problem.

A1,A2,A3,A4,A6 are fixed, three dimensional coordinates associated with each position on the ring.

Bonds = str(input('Enter pi bond atomic pair in format "A1-A2" or type PLOT: '))

while (Bonds != 'PLOT'):
    Bonds2 = Bonds.split('-')
    for bond in Bonds2:
        if bond == 'A1':
            pi_list.append(A1)
        elif bond == 'A2':
            pi_list.append(A2)
        elif bond == 'A3':
            pi_list.append(A3)
        elif bond == 'A4':
            pi_list.append(A4)
        elif bond == 'A5':
            pi_list.append(A5)
        elif bond == 'A6':
            pi_list.append(A6)
    if pi_list != []:
        x += pi_bond(x,pi_list[0], pi_list[1])
    for atom in Coords:
        for pi_atom in pi_list:
        if (atom != pi_list[0]) and (atom != pi_list[1]) and (float(distance(pi_atom,atom))==float(distance(A1,A2))):
            x += single_bond(pi_atom, atom)
    pi_list.clear()

    if (Bonds == 'PLOT'):
        break
    Bonds = str(input('Enter additional pi bonds or type PLOT: '))

Can anyone think of a more elegant way to approach this?

Here is the function I am currently using to add a double bond between two points:

def pi_bond(x,Cx,Cy):
    X = 0.05
    if Cx[1] == Cy[1]:
        Ca = Cx[1]-X
        Cb = Cx[1]+X
        Cc = Cy[1]-X
        Cd = Cy[1]+X
        Cx1 = (Cx[0], Ca, Cx[2])
        Cx2 = (Cx[0], Cb, Cx[2])
        Cy1 = (Cy[0], Cc, Cy[2])
        Cy2 = (Cy[0], Cd, Cy[2])
        x += LineSegment(Cx1, Cy1, 1, color='white', axes=False, frame=False)
        x += LineSegment(Cx2, Cy2, 1, color='white', axes=False, frame=False)
       return x
    else:
        Ca = Cx[2]-X
        Cb = Cx[2]+X
        Cc = Cy[2]-X
        Cd = Cy[2]+X
        Cx1 = (Cx[0], Cx[1], Ca)
        Cx2 = (Cx[0], Cx[1], Cb)
        Cy1 = (Cy[0], Cy[1], Cc)
        Cy2 = (Cy[0], Cy[1], Cd)
        x += LineSegment(Cx1, Cy1, 1, color='white', axes=False, frame=False)
        x += LineSegment(Cx2, Cy2, 1, color='white', axes=False, frame=False)
        return x
2020-03-14 12:52:36 -0600 commented answer How to import implicit_plot3d into .py file?

Yes, Normalize adjusts all the radii relative to a fixed size for carbon. Thank you. This is exactly what I was looking for.

2020-03-14 12:43:15 -0600 received badge  Editor (source)
2020-03-14 12:28:26 -0600 asked a question How to import implicit_plot3d into .py file?

I'm working on some code for creating 3D molecular structures and everything seems to be working well aside from the fact that I have no idea how to import the 'implicit_plot3d' function into my script. I want to run the code as 'Molecule.py' directly inside of sage. When I try to run it right now I get the following error:

NameError: name 'implicit_plot3d' is not defined

Can anyone help me identify the issue? I'm not very experienced with either Python or Sage.

import sys
from sage.plot.plot3d.shapes import LineSegment, Sphere
from sage.all import rainbow

BC = ['white','black','blue','red','grey','red','pink','cyan']

Radii = {'O': 0.73, 'N': 0.75, 'C': 0.77, 'He': 0.32, 'H': 0.37, 'S': 1.02, 'Cl': 0.99, 
'F': 0.71, 'Xe': 1.30, 'Si': 1.11, 'B': 0.82, 'P': 1.06, 'Br': 1.14}

Valency = {'O': 6, 'N': 5, 'C': 4, 'He': 0, 'H': 1, 'S': 6, 'Cl': 7, 
'F': 7, 'Xe': 8, 'Si': 4, 'B': 3, 'P': 5, 'Br': 7}

Colors = {'O': 'red', 'N': 'green', 'C': 'black', 'He': 'cyan', 'H': 'white', 'S': 'yellow',
'Si': 'purple', 'Xe': 'pink', 'F': 'blue', 'Cl': 'green', 'B': 'pink', 'P': 'orange',
'Br': 'red'}

Here's a sample function

def bent_geometry():

    x = input("Enter bent molecule A: ")
    y = input("Enter bent molecule B: ")
    X =(0,0,0)
    Y1 = (0, 0.75, -0.70)
    Y2 = (0, -0.75,-0.70)

    C2_a = (0, 0, -1.5)
    C2_b = (0, 0, 1.5)

    XY2  = Sphere(Normalize(Radii[x]), color=Colors[x]).translate(X)
    XY2 += Sphere(Normalize(Radii[y]), color=Colors[y]).translate(Y1)
    XY2 += Sphere(Normalize(Radii[y]), color=Colors[y]).translate(Y2)

    XY2 += LineSegment(X, Y1, 1, color=BC[0], axes=False, frame=False)
    XY2 += LineSegment(X, Y2, 1, color=BC[0], axes=False, frame=False)
    XY2 += LineSegment(C2_a, C2_b, 1, color = BC[2], axes=True, frame=True)

    XY2 += implicit_plot3d(lambda x,y,z: x, (-1,1), (-1,1), (-1,1),color='cyan',opacity=0.8)
    XY2 += implicit_plot3d(lambda x,y,z: y, (-1,1), (-1,1), (-1,1),color='white',opacity=0.8)

    XY2 += text3d('C2',(0,0,1.575))
    return XY2
2020-03-08 11:37:12 -0600 commented answer How to obtain decimal value from .solve() ?

Thank for your such a thorough explanation! As someone new to Sage it means a lot when others are willing to break down and explain trivial issues like this.

2020-03-07 09:10:35 -0600 asked a question How to obtain decimal value from .solve() ?
var('n')
f(x,n) = -x^-1 + (1.097e7)*((1/(n-1)^2) - (1/n^2))

for i in range (2,10):
    L = f(x,i).solve(x)
    print(L)

[x == (1/8227500)]

[x == (9/13712500)]

....

How can I have this display a float value for x, rather than a fraction?