Changeset 1948
- Timestamp:
- Aug 3, 2015 8:33:08 AM (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIspc.py
r1917 r1948 1472 1472 return A 1473 1473 1474 def DoFrac(): 1475 delt2 = np.eye(2)*0.001 1476 FSC = np.ones(2,dtype='i') 1477 VFSC = np.ones(2) 1478 CSI = [np.zeros((2),dtype='i'),np.zeros(2)] 1479 if 'Crenel' in waveType: 1480 dF = fracCrenel(tau,delt2[:1],delt2[1:]).squeeze() 1481 else: 1482 dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze() 1483 dFT = np.zeros_like(dF) 1484 dFTP = [] 1485 for i in SdIndx: 1486 sop = Sop[i] 1487 ssop = SSop[i] 1488 ssopinv = nl.inv(ssop[0]) 1489 mst = ssopinv[3][:3] 1490 epsinv = ssopinv[3][3] 1491 sdet = nl.det(sop[0]) 1492 ssdet = nl.det(ssop[0]) 1493 dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3] 1494 dT = 1.0 1495 if np.any(dtau%.5): 1496 dT = np.tan(np.pi*np.sum(dtau%.5)) 1497 tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]) 1498 fsc = np.ones(2,dtype='i') 1499 if 'Crenel' in waveType: 1500 dFT = fracCrenel(tauT,delt2[:1],delt2[1:]).squeeze() 1501 fsc = [1,1] 1502 else: #Fourier 1503 dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze() 1504 dFT = nl.det(sop[0])*dFT 1505 dFT = dFT[:,np.argsort(tauT)] 1506 dFT[0] *= ssdet 1507 dFT[1] *= sdet 1508 dFTP.append(dFT) 1509 1510 if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']): 1511 fsc = [1,1] 1512 CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]] 1513 FSC = np.zeros(2,dtype='i') 1514 return CSI,dF,dFTP 1515 else: 1516 for i in range(2): 1517 if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6): 1518 fsc[i] = 1 1519 else: 1520 fsc[i] = 0 1521 FSC &= fsc 1522 print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc 1523 n = -1 1524 for i,F in enumerate(FSC): 1525 if F: 1526 n += 1 1527 CSI[0][i] = n+1 1528 CSI[1][i] = 1.0 1529 1530 return CSI,dF,dFTP 1531 1532 def DoXYZ(): 1533 delt6 = np.eye(6)*0.001 1534 if 'Fourier' in waveType: 1535 dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis] 1536 #3x6x12 modulated position array (X,Spos,tau)& force positive 1537 CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))] 1538 elif waveType == 'Sawtooth': 1539 CSI = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])] 1540 elif waveType == 'ZigZag': 1541 CSI = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])] 1542 XSC = np.ones(6,dtype='i') 1543 dXTP = [] 1544 for i in SdIndx: 1545 sop = Sop[i] 1546 ssop = SSop[i] 1547 xsc = np.ones(6,dtype='i') 1548 ssopinv = nl.inv(ssop[0]) 1549 mst = ssopinv[3][:3] 1550 epsinv = ssopinv[3][3] 1551 sdet = nl.det(sop[0]) 1552 ssdet = nl.det(ssop[0]) 1553 dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3] 1554 dT = 1.0 1555 if np.any(dtau%.5): 1556 dT = np.tan(np.pi*np.sum(dtau%.5)) 1557 tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]) 1558 if waveType == 'Fourier': 1559 dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis] 1560 elif waveType == 'Sawtooth': 1561 dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis] 1562 elif waveType == 'ZigZag': 1563 dXT = posZigZag(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis] 1564 dXT = np.inner(sop[0],dXT.T) 1565 dXT = np.swapaxes(dXT,1,2) 1566 dXT[:,:3,:] *= ssdet 1567 dXTP.append(dXT) 1568 if waveType == 'Fourier': 1569 if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']): 1570 xsc[3:6] = 0 1571 CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]], 1572 [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]] 1573 if '(x)' in siteSym: 1574 CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.] 1575 if 'm' in siteSym and len(SdIndx) == 1: 1576 CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.] 1577 elif '(y)' in siteSym: 1578 CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.] 1579 if 'm' in siteSym and len(SdIndx) == 1: 1580 CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.] 1581 elif '(z)' in siteSym: 1582 CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.] 1583 if 'm' in siteSym and len(SdIndx) == 1: 1584 CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.] 1585 elif '(xy)' in siteSym: 1586 CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]] 1587 CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]] 1588 elif np.allclose(dX[0,0,:],dXT[0,1,:]*sdet): 1589 if 'xy' in siteSym or '+-0' in siteSym: 1590 CSI[0][3:5] = [[12,0,0],[12,0,0]] 1591 CSI[1][3:5] = [[1.,0,0],[-sdet,0,0]] 1592 xsc[3:5] = 0 1593 else: 1594 for i in range(3): 1595 if np.allclose(dX[i,i,:],dXT[i,i,:]*sdet): 1596 xsc[i] = 1 1597 else: 1598 xsc[i] = 0 1599 if np.allclose(dX[i,i+3,:],dXT[i,i+3,:]): 1600 xsc[i+3] = 1 1601 else: 1602 xsc[i+3] = 0 1603 XSC &= xsc 1604 print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc 1605 n = -1 1606 print XSC 1607 for i,X in enumerate(XSC): 1608 if X: 1609 n += 1 1610 CSI[0][i][0] = n+1 1611 CSI[1][i][0] = 1.0 1612 1613 return CSI,dX,dXTP 1614 1615 def DoUij(): 1616 delt12 = np.eye(12)*0.0001 1617 dU = posFourier(tau,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x12 array 1618 CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))] 1619 USC = np.ones(12,dtype='i') 1620 dUTP = [] 1621 for i in SdIndx: 1622 sop = Sop[i] 1623 ssop = SSop[i] 1624 ssopinv = nl.inv(ssop[0]) 1625 mst = ssopinv[3][:3] 1626 epsinv = ssopinv[3][3] 1627 sdet = nl.det(sop[0]) 1628 ssdet = nl.det(ssop[0]) 1629 dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3] 1630 dT = 1.0 1631 if np.any(dtau%.5): 1632 dT = np.tan(np.pi*np.sum(dtau%.5)) 1633 tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]) 1634 usc = np.ones(12,dtype='i') 1635 dUT = posFourier(tauT,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x49 array 1636 dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3) #convert dUT to 12x49x3x3 1637 dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) 1638 dUT = np.array(U2Uij(dUijT)) 1639 dUT = dUT[:,:,np.argsort(tauT)] 1640 dUTP.append(dUT) 1641 if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']): 1642 CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 1643 [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]], 1644 [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.], 1645 [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]] 1646 if 'mm2(x)' in siteSym: 1647 CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.] 1648 USC = [1,1,1,0,1,0,1,1,1,0,1,0] 1649 elif '(xy)' in siteSym: 1650 CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0], 1651 [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]] 1652 CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]] 1653 USC = [1,1,1,1,1,1,1,1,1,1,1,1] 1654 elif '(x)' in siteSym: 1655 CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.] 1656 elif '(y)' in siteSym: 1657 CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.] 1658 elif '(z)' in siteSym: 1659 CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.] 1660 for i in range(6): 1661 if not USC[i]: 1662 CSI[0][i] = [0,0,0] 1663 CSI[1][i] = [0.,0.,0.] 1664 CSI[0][i+6] = [0,0,0] 1665 CSI[1][i+6] = [0.,0.,0.] 1666 else: 1667 for i in range(6): 1668 if np.allclose(dU[i,i,:],-dUT[i,i,:]): 1669 usc[i] = 1 1670 else: 1671 usc[i] = 0 1672 if np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): 1673 usc[i+6] = 1 1674 else: 1675 usc[i+6] = 0 1676 if np.any(dUT[0,1,:]): 1677 if '4/m' in siteSym: 1678 CSI[0][6:8] = [[12,0,0],[12,0,0]] 1679 if ssop[1][3]: 1680 CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]] 1681 usc[9] = 1 1682 else: 1683 CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]] 1684 usc[9] = 0 1685 elif '4' in siteSym: 1686 CSI[0][6:8] = [[12,0,0],[12,0,0]] 1687 CSI[0][:2] = [[11,0,0],[11,0,0]] 1688 if ssop[1][3]: 1689 CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]] 1690 CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]] 1691 usc[2] = 0 1692 usc[8] = 0 1693 usc[3] = 1 1694 usc[9] = 1 1695 else: 1696 CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]] 1697 CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]] 1698 usc[2] = 1 1699 usc[8] = 1 1700 usc[3] = 0 1701 usc[9] = 0 1702 elif 'xy' in siteSym or '+-0' in siteSym: 1703 if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet): 1704 print np.allclose(dU[0,0,:],dUT[0,1,:]*sdet),sdet 1705 CSI[0][4:6] = [[12,0,0],[12,0,0]] 1706 CSI[0][6:8] = [[11,0,0],[11,0,0]] 1707 CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]] 1708 CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]] 1709 usc[4:6] = 0 1710 usc[6:8] = 0 1711 1712 print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc 1713 USC &= usc 1714 if not np.any(dtau%.5): 1715 n = -1 1716 for i,U in enumerate(USC): 1717 if U: 1718 n += 1 1719 CSI[0][i][0] = n+1 1720 CSI[1][i][0] = 1.0 1721 1722 return CSI,dU,dUTP 1723 1724 def DoUiso(): 1725 print 'Uiso' 1726 delt4 = np.eye(4)*0.001 1727 1474 1728 print 'super space group: ',SSGData['SSpGrp'] 1475 1729 CSI = {'Sfrac':[[[1,0],[2,0]],[[1.,0.],[1.,0.]]], … … 1507 1761 SGOps.append([-op[0],-op[1]%1.]) 1508 1762 SSGOps.append([-sop[0],-sop[1]%1.]) 1509 #build set of sym ops around special po asition1763 #build set of sym ops around special position 1510 1764 SSop = [] 1511 1765 Sop = [] … … 1525 1779 print 'special pos super operators: ',SSOpText 1526 1780 #setup displacement arrays 1527 tau = np.linspace(0,1,49,True) 1528 delt2 = np.eye(2)*0.001 1529 delt4 = np.eye(4)*0.001 1530 delt6 = np.eye(6)*0.001 1531 delt12 = np.eye(12)*0.0001 1781 tau = np.linspace(-1,1,49,True) 1532 1782 #make modulation arrays - one parameter at a time 1533 1783 #site fractions 1534 CSI['Sfrac'] = [np.zeros((2),dtype='i'),np.ones(2)] 1535 if 'Crenel' in waveType: 1536 dF = fracCrenel(tau,delt2[:1],delt2[1:]).squeeze() 1537 else: 1538 dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze() 1539 dFT = np.zeros_like(dF) 1540 #positions 1541 if 'Fourier' in waveType: 1542 dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis] 1543 #3x6x12 modulated position array (X,Spos,tau)& force positive 1544 CSI['Spos'] = [np.zeros((6,3),dtype='i'),np.zeros((6,3))] 1545 elif waveType == 'Sawtooth': 1546 CSI['Spos'] = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])] 1547 elif waveType == 'ZigZag': 1548 CSI['Spos'] = [np.array([[1,],[2,],[3,],[4,]]),np.array([[1.0,],[1.0,],[1.0,],[1.0,]])] 1784 CSI['Sfrac'],dF,dFTP = DoFrac() 1785 #positions 1786 CSI['Spos'],dX,dXTP = DoXYZ() 1549 1787 #anisotropic thermal motion 1550 dU = posFourier(tau,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x12 array 1551 CSI['Sadp'] = [np.zeros((12,3),dtype='i'),np.zeros((12,3))] 1552 1553 FSC = np.ones(2,dtype='i') 1554 VFSC = np.ones(2) 1555 XSC = np.ones(6,dtype='i') 1556 USC = np.ones(12,dtype='i') 1557 dFTP = [] 1558 dXTP = [] 1559 dUTP = [] 1560 for i in SdIndx: 1561 sop = Sop[i] 1562 ssop = SSop[i] 1563 fsc = np.ones(2,dtype='i') 1564 xsc = np.ones(6,dtype='i') 1565 ssopinv = nl.inv(ssop[0]) 1566 mst = ssopinv[3][:3] 1567 epsinv = ssopinv[3][3] 1568 sdet = nl.det(sop[0]) 1569 ssdet = nl.det(ssop[0]) 1570 dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3] 1571 dT = 1.0 1572 if np.any(dtau%.5): 1573 dT = np.tan(np.pi*np.sum(dtau%.5)) 1574 tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]) 1575 if waveType == 'Fourier': 1576 dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis] 1577 elif waveType == 'Sawtooth': 1578 dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis] 1579 elif waveType == 'ZigZag': 1580 dXT = posZigZag(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis] 1581 dXT = np.inner(sop[0],dXT.T) 1582 dXT = np.swapaxes(dXT,1,2) 1583 dXT[:,:3,:] *= ssdet 1584 dXTP.append(dXT) 1585 if waveType == 'Fourier': 1586 if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']): 1587 CSI['Spos'] = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]], 1588 [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]] 1589 if '(x)' in siteSym: 1590 CSI['Spos'][1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.] 1591 if 'm' in siteSym and len(SdIndx) == 1: 1592 CSI['Spos'][1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.] 1593 elif '(y)' in siteSym: 1594 CSI['Spos'][1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.] 1595 if 'm' in siteSym and len(SdIndx) == 1: 1596 CSI['Spos'][1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.] 1597 elif '(z)' in siteSym: 1598 CSI['Spos'][1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.] 1599 if 'm' in siteSym and len(SdIndx) == 1: 1600 CSI['Spos'][1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.] 1601 elif '(xy)' in siteSym: 1602 CSI['Spos'][0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]] 1603 CSI['Spos'][1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]] 1604 for i in range(3): 1605 if not XSC[i]: 1606 CSI['Spos'][0][i] = [0,0,0] 1607 CSI['Spos'][1][i] = [0.,0.,0.] 1608 CSI['Spos'][0][i+3] = [0,0,0] 1609 CSI['Spos'][1][i+3] = [0.,0.,0.] 1610 elif np.allclose(dX[0,0,:],dXT[0,1,:]*sdet): 1611 if 'xy' in siteSym or '+-0' in siteSym: 1612 CSI['Spos'][0][3:5] = [[12,0,0],[12,0,0]] 1613 CSI['Spos'][1][3:5] = [[1.,0,0],[-sdet,0,0]] 1614 xsc[3:5] = 0 1615 else: 1616 for i in range(3): 1617 if np.allclose(dX[i,i,:],dXT[i,i,:]*sdet): 1618 xsc[i] = 1 1619 else: 1620 xsc[i] = 0 1621 if np.allclose(dX[i,i+3,:],dXT[i,i+3,:]): 1622 xsc[i+3] = 1 1623 else: 1624 xsc[i+3] = 0 1625 XSC &= xsc 1626 1627 fsc = np.ones(2,dtype='i') 1628 if 'Crenel' in waveType: 1629 dFT = fracCrenel(tauT,delt2[:1],delt2[1:]).squeeze() 1630 fsc = [1,1] 1631 else: 1632 dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze() 1633 dFT = nl.det(sop[0])*dFT 1634 dFT = dFT[:,np.argsort(tauT)] 1635 dFT[0] *= ssdet 1636 dFT[1] *= sdet 1637 dFTP.append(dFT) 1638 1639 if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']): 1640 fsc = [1,1] 1641 CSI['Sfrac'] = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]] 1642 for i in range(2): 1643 if not FSC[i]: 1644 CSI['Sfrac'][0][i] = [0,0] 1645 CSI['Sfrac'][1][i] = [0.,0.] 1646 else: 1647 for i in range(2): 1648 if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6): 1649 fsc[i] = 1 1650 else: 1651 fsc[i] = 0 1652 FSC &= fsc 1653 1654 usc = np.ones(12,dtype='i') 1655 dUT = posFourier(tauT,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x49 array 1656 dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3) #convert dUT to 12x49x3x3 1657 dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) 1658 dUT = np.array(U2Uij(dUijT)) 1659 dUT = dUT[:,:,np.argsort(tauT)] 1660 dUTP.append(dUT) 1661 if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']): 1662 CSI['Sadp'] = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0], 1663 [1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]], 1664 [[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.], 1665 [1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]] 1666 if 'mm2(x)' in siteSym: 1667 CSI['Sadp'][1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.] 1668 USC = [1,1,1,0,1,0,1,1,1,0,1,0] 1669 elif '(xy)' in siteSym: 1670 CSI['Sadp'][0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0], 1671 [1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]] 1672 CSI['Sadp'][1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]] 1673 USC = [1,1,1,1,1,1,1,1,1,1,1,1] 1674 elif '(x)' in siteSym: 1675 CSI['Sadp'][1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.] 1676 elif '(y)' in siteSym: 1677 CSI['Sadp'][1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.] 1678 elif '(z)' in siteSym: 1679 CSI['Sadp'][1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.] 1680 for i in range(6): 1681 if not USC[i]: 1682 CSI['Sadp'][0][i] = [0,0,0] 1683 CSI['Sadp'][1][i] = [0.,0.,0.] 1684 CSI['Sadp'][0][i+6] = [0,0,0] 1685 CSI['Sadp'][1][i+6] = [0.,0.,0.] 1686 else: 1687 for i in range(6): 1688 if np.allclose(dU[i,i,:],-dUT[i,i,:]): 1689 usc[i] = 1 1690 else: 1691 usc[i] = 0 1692 if np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): 1693 usc[i+6] = 1 1694 else: 1695 usc[i+6] = 0 1696 if np.any(dUT[0,1,:]): 1697 if '4/m' in siteSym: 1698 CSI['Sadp'][0][6:8] = [[12,0,0],[12,0,0]] 1699 if ssop[1][3]: 1700 CSI['Sadp'][1][6:8] = [[1.,0.,0.],[-1.,0.,0.]] 1701 usc[9] = 1 1702 else: 1703 CSI['Sadp'][1][6:8] = [[1.,0.,0.],[1.,0.,0.]] 1704 usc[9] = 0 1705 elif '4' in siteSym: 1706 CSI['Sadp'][0][6:8] = [[12,0,0],[12,0,0]] 1707 CSI['Sadp'][0][:2] = [[11,0,0],[11,0,0]] 1708 if ssop[1][3]: 1709 CSI['Sadp'][1][:2] = [[1.,0.,0.],[-1.,0.,0.]] 1710 CSI['Sadp'][1][6:8] = [[1.,0.,0.],[-1.,0.,0.]] 1711 usc[2] = 0 1712 usc[8] = 0 1713 usc[3] = 1 1714 usc[9] = 1 1715 else: 1716 CSI['Sadp'][1][:2] = [[1.,0.,0.],[1.,0.,0.]] 1717 CSI['Sadp'][1][6:8] = [[1.,0.,0.],[1.,0.,0.]] 1718 usc[2] = 1 1719 usc[8] = 1 1720 usc[3] = 0 1721 usc[9] = 0 1722 elif 'xy' in siteSym or '+-0' in siteSym: 1723 if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet): 1724 print np.allclose(dU[0,0,:],dUT[0,1,:]*sdet),sdet 1725 CSI['Sadp'][0][4:6] = [[12,0,0],[12,0,0]] 1726 CSI['Sadp'][0][6:8] = [[11,0,0],[11,0,0]] 1727 CSI['Sadp'][1][4:6] = [[1.,0.,0.],[sdet,0.,0.]] 1728 CSI['Sadp'][1][6:8] = [[1.,0.,0.],[sdet,0.,0.]] 1729 usc[4:6] = 0 1730 usc[6:8] = 0 1731 1732 print SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc 1733 USC &= usc 1734 if not np.any(dtau%.5): 1735 n = -1 1736 for i,U in enumerate(USC): 1737 if U: 1738 n += 1 1739 CSI['Sadp'][0][i][0] = n+1 1740 CSI['Sadp'][1][i][0] = 1.0 1741 if waveType == 'Fourier': 1742 n = -1 1743 for i,X in enumerate(XSC): 1744 if X: 1745 n += 1 1746 CSI['Spos'][0][i][0] = n+1 1747 CSI['Spos'][1][i][0] = 1.0 1748 n = -1 1749 for i,F in enumerate(FSC): 1750 if F: 1751 n += 1 1752 CSI['Sfrac'][0][i] = n+1 1753 CSI['Sfrac'][1][i] = 1.0 1754 else: 1755 CSI['Sfrac'][0][i] = 0 1756 CSI['Sfrac'][1][i] = 0. 1788 CSI['Sadp'],dU,dUTP = DoUij() 1757 1789 CSI['Spos'][0] = orderParms(CSI['Spos'][0]) 1758 1790 CSI['Sadp'][0] = orderParms(CSI['Sadp'][0]) 1759 1791 if debug: 1760 return CSI,[tau,tau T],[dF,dFTP],[dX,dXTP],[dU,dUTP]1792 return CSI,[tau,tau],[dF,dFTP],[dX,dXTP],[dU,dUTP] 1761 1793 else: 1762 1794 return CSI
Note: See TracChangeset
for help on using the changeset viewer.