Changeset 1927
- Timestamp:
- Jul 14, 2015 3:55:03 PM (8 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r1926 r1927 376 376 return XYZ 377 377 378 def FindNeighbors(phase,FrstName,AtNames ):378 def FindNeighbors(phase,FrstName,AtNames,notName=''): 379 379 General = phase['General'] 380 380 cx,ct,cs,cia = General['AtomPtrs'] … … 400 400 for j in IndB[0]: 401 401 if j != Orig: 402 Neigh.append([AtNames[j],dist[j]]) 403 Ids.append(Atoms[j][cia+8]) 402 if AtNames[j] != notName: 403 Neigh.append([AtNames[j],dist[j]]) 404 Ids.append(Atoms[j][cia+8]) 404 405 return Neigh,[OId,Ids] 405 406 406 407 def AddHydrogens(AtLookUp,General,Atoms,AddHydId): 408 407 409 cx,ct,cs,cia = General['AtomPtrs'] 408 410 Cell = General['Cell'][1:7] 409 411 Amat,Bmat = G2lat.cell2AB(Cell) 410 nBonds = AddHydId[ 2]+len(AddHydId[1])412 nBonds = AddHydId[-1]+len(AddHydId[1]) 411 413 Oatom = GetAtomsById(Atoms,AtLookUp,[AddHydId[0],])[0] 412 414 OXYZ = np.array(Oatom[cx:cx+3]) … … 415 417 DX = np.inner(Amat,TXYZ-OXYZ).T 416 418 if nBonds == 4: 417 if AddHydId[ 2] == 1:419 if AddHydId[-1] == 1: 418 420 Vec = TXYZ-OXYZ 419 421 Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2,axis=0)) … … 422 424 Hpos = OXYZ-0.98*np.inner(Bmat,Vec).T/Len 423 425 return [Hpos,] 424 elif AddHydId[2] == 2: 425 print 'add 2 H' 426 return [] 426 elif AddHydId[-1] == 2: 427 Vec = np.inner(Amat,TXYZ-OXYZ).T 428 Vec[0] += Vec[1] #U - along bisector 429 Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,np.newaxis] 430 Mat2 = np.cross(Vec[0],Vec[1]) #UxV 431 Mat2 /= np.sqrt(np.sum(Mat2**2)) 432 Mat3 = np.cross(Mat2,Vec[0]) #(UxV)xU 433 iMat = nl.inv(np.array([Vec[0],Mat2,Mat3])) 434 Hpos = np.array([[-0.97*cosd(54.75),0.97*sind(54.75),0.], 435 [-0.97*cosd(54.75),-0.97*sind(54.75),0.]]) 436 Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ 437 return Hpos 427 438 else: 428 print 'add 3 H' 429 return [] 439 Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0] 440 RXYZ = np.array(Ratom[cx:cx+3]) 441 Vec = np.inner(Amat,np.array([OXYZ-TXYZ[0],RXYZ-TXYZ[0]])).T 442 Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,np.newaxis] 443 Mat2 = np.cross(Vec[0],Vec[1]) #UxV 444 Mat2 /= np.sqrt(np.sum(Mat2**2)) 445 Mat3 = np.cross(Mat2,Vec[0]) #(UxV)xU 446 iMat = nl.inv(np.array([Vec[0],Mat2,Mat3])) 447 a = 0.96*cosd(70.5) 448 b = 0.96*sind(70.5) 449 Hpos = np.array([[a,0.,-b],[a,-b*cosd(30.),0.5*b],[a,b*cosd(30.),0.5*b]]) 450 Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ 451 return Hpos 430 452 elif nBonds == 3: 431 if AddHydId[ 2] == 1:453 if AddHydId[-1] == 1: 432 454 Vec = np.sum(TXYZ-OXYZ,axis=0) 433 455 Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2)) … … 435 457 Hpos = OXYZ+Vec 436 458 return [Hpos,] 437 elif AddHydId[2] == 2: 459 elif AddHydId[-1] == 2: 460 Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0] 461 RXYZ = np.array(Ratom[cx:cx+3]) 462 Vec = np.inner(Amat,np.array([OXYZ-TXYZ[0],RXYZ-TXYZ[0]])).T 463 Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,np.newaxis] 464 Mat2 = np.cross(Vec[0],Vec[1]) #UxV 465 Mat2 /= np.sqrt(np.sum(Mat2**2)) 466 Mat3 = np.cross(Mat2,Vec[0]) #(UxV)xU 467 iMat = nl.inv(np.array([Vec[0],Mat2,Mat3])) 438 468 print 'add 2 H' 439 469 return [] 440 470 else: #2 bonds 441 471 if 'C' in Oatom[ct]: 442 print 'sp atom',Oatom[ct-1] 443 return [] 472 Vec = TXYZ[0]-OXYZ 473 Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2)) 474 Vec = -0.93*Vec/Len 475 Hpos = OXYZ+Vec 476 return [Hpos,] 444 477 elif 'O' in Oatom[ct]: 445 #idea - construct ring at 0.82 from O atom & find high spot? 446 print 'sp3 atom ',Oatom[ct-1] 447 print 'add 1 H' 448 return [] 478 mapData = General['Map'] 479 Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0] 480 RXYZ = np.array(Ratom[cx:cx+3]) 481 Vec = np.inner(Amat,np.array([OXYZ-TXYZ[0],RXYZ-TXYZ[0]])).T 482 Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,np.newaxis] 483 Mat2 = np.cross(Vec[0],Vec[1]) #UxV 484 Mat2 /= np.sqrt(np.sum(Mat2**2)) 485 Mat3 = np.cross(Mat2,Vec[0]) #(UxV)xU 486 iMat = nl.inv(np.array([Vec[0],Mat2,Mat3])) 487 a = 0.82*cosd(70.5) 488 b = 0.82*sind(70.5) 489 azm = np.arange(0.,360.,5.) 490 Hpos = np.array([[a,b*cosd(x),b*sind(x)] for x in azm]) 491 Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ 492 return Hpos 449 493 return [] 450 451 452 453 454 494 455 495 def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj): #unfinished & not used -
trunk/GSASIIphsGUI.py
r1926 r1927 1572 1572 nH = 4-len(neigh[1][0]) 1573 1573 bonds = dict(neigh[1][0]) 1574 nextName = '' 1575 if len(bonds) == 1: 1576 nextName = bonds.keys()[0] 1574 1577 for bond in bonds: 1575 1578 if 'C' in neigh[0]: … … 1583 1586 nH -= 1 1584 1587 break 1588 nextneigh = [] 1589 if nextName: 1590 nextneigh = G2mth.FindNeighbors(data,nextName,AtNames,notName=neigh[0]) 1591 neigh[1][1].append(nextneigh[1][1][0]) 1585 1592 neigh[2] = max(0,nH) #set expected no. H's needed 1586 1593 AddHydIds.append(neigh[1][1]) … … 1591 1598 Nat = len(atomData) 1592 1599 Neigh = dlg.GetData() 1593 for ineigh,neigh in enumerate(Neigh): 1600 mapData = generalData['Map'] 1601 mapError = False 1602 for ineigh,neigh in enumerate(Neigh): 1594 1603 AddHydIds[ineigh].append(neigh[2]) 1604 if 'O' in neigh[0] and not len(mapData['rho']) or not 'delt-F' in mapData['MapType']: 1605 mapError = True 1606 continue 1595 1607 Hxyz = G2mth.AddHydrogens(AtLookUp,generalData,atomData,AddHydIds[ineigh]) 1596 1608 for X in Hxyz: 1597 1609 Nat += 1 1598 1610 AtomAdd(X[0],X[1],X[2],'H','H(%d)'%(Nat)) 1599 1611 if mapError: 1612 G2frame.ErrorDialog('Add H atom error','Adding O-H atoms requires delt-F map') 1600 1613 SetupGeneral() 1601 1614 FillAtomsGrid(Atoms)
Note: See TracChangeset
for help on using the changeset viewer.