Changeset 965
 Timestamp:
 Jun 24, 2013 12:20:48 PM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r962 r965 32 32 import GSASIIspc as G2spc 33 33 import numpy.fft as fft 34 import pypowder as pyd 34 35 35 36 sind = lambda x: np.sin(x*np.pi/180.) … … 2317 2318 x0 = schedule.getstart_temp(best_state) 2318 2319 else: 2319 x0 = random.uniform(size=len(x0))*(upperlower) + lower2320 # x0 = random.uniform(size=len(x0))*(upperlower) + lower 2320 2321 best_state.x = None 2321 2322 best_state.cost = numpy.Inf … … 2506 2507 Mdata.append(float(len(G2spc.GenAtom(xyz,SGData)))) 2507 2508 return np.array(Mdata) 2508 2509 2509 2510 def GetAtomTX(RBdata,parmDict): 2510 2511 'Needs a doc string' … … 2519 2520 parm = ':'+str(iatm)+key 2520 2521 if parm in parmDict: 2521 keys[key][iatm] = parmDict[parm] 2522 if key == ':Atype': 2523 keys[key][iatm] = aTypes.index(parmDict[parm]) 2524 else: 2525 keys[key][iatm] = parmDict[parm] 2522 2526 iatm = nfixAt 2523 2527 for iObj in range(parmDict['nObj']): … … 2526 2530 if parmDict[pfx+'Type'] == 'Vector': 2527 2531 RBRes = RBdata['Vector'][parmDict[pfx+'RBId']] 2528 aTypes = RBRes['rbTypes']2529 2532 vecs = RBRes['rbVect'] 2530 2533 mags = RBRes['VectMag'] … … 2534 2537 elif parmDict[pfx+'Type'] == 'Residue': 2535 2538 RBRes = RBdata['Residue'][parmDict[pfx+'RBId']] 2536 aTypes = RBRes['rbTypes']2537 2539 Cart = np.array(RBRes['rbXYZ']) 2538 2540 for itor,seq in enumerate(RBRes['rbSeq']): … … 2549 2551 for j in range(3): 2550 2552 Xdata[j][iatm] = X[j] 2551 Tdata[iatm] = aTypes [i]2553 Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i]) 2552 2554 iatm += 1 2553 2555 elif parmDict[pfx+'Type'] == 'Atom': … … 2556 2558 parm = pfx+key[1:] #remove extra ':' 2557 2559 if parm in parmDict: 2558 keys[key][atNo] = parmDict[parm] 2560 if key == ':Atype': 2561 keys[key][atNo] = aTypes.index(parmDict[parm]) 2562 else: 2563 keys[key][atNo] = parmDict[parm] 2559 2564 iatm += 1 2560 2565 else: … … 2563 2568 2564 2569 def calcMDcorr(MDval,MDaxis,Uniq,G): 2565 sumMD = 0 2566 for H in Uniq: 2567 cosP,sinP = G2lat.CosSinAngle(H,MDaxis,G) 2568 A = 1.0/np.sqrt((MDval*cosP)**2+sinP**2/MDval) 2569 sumMD += A**3 2570 return sumMD 2570 sumMD = len(Uniq) 2571 if MDval != 1.0: 2572 for H in Uniq: 2573 cosP,sinP = G2lat.CosSinAngle(H,MDaxis,G) 2574 A = 1.0/np.sqrt((MDval*cosP)**2+sinP**2/MDval) 2575 sumMD += A**3 2576 sumMD = np.sum(1./np.sqrt((MDval*cosP)**2+sinP**2/MDval)**3) 2577 return sumMD/len(Uniq) 2578 2579 def mcsasfCalc(ifInv,Tdata,Mdata,Xdata,mul,FFs,Uniq,Phi): 2580 ''' Code here needs to be in fortran''' 2581 Icalc = pyd.pymcsasfcalc(ifInv,len(Tdata),Tdata,Mdata,Xdata.flatten(), 2582 mul,len(FFs),FFs,len(Uniq),Uniq.flatten(),Phi) 2583 return Icalc 2584 # FF = np.zeros(len(Tdata)) 2585 # for i,j in enumerate(Tdata): #NB: generator here is slower! 2586 # FF[i] = FFs[j] 2587 # FF *= Mdata/len(Phi) #FF*occ 2588 # phase = np.inner(Uniq,Xdata) #hx+ky+lz 2589 # phase += Phi[:,np.newaxis] #hx+ky+lz+phi 2590 # if len(Uniq) == 3: print 'python',phase 2591 # cosp = np.cos(twopi*phase) 2592 # fas = np.sum(FF*cosp) 2593 # print fas 2594 # if ifInv: 2595 # fbs = 0. 2596 # else: 2597 # sinp = np.sin(twopi*phase) 2598 # fbs = np.sum(FF*sinp) 2599 # return (fas**2+fbs**2)*mul 2571 2600 2572 2601 def mcsaCalc(values,refList,rcov,ifInv,RBdata,varyList,parmDict): 2573 #''' Compute structure factors for all h,k,l for phase2574 #input:2575 #refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]2576 #ParmDict:2577 #puts result F^2 in each ref[8] in refList2578 #'''2602 ''' Compute structure factors for all h,k,l for phase 2603 input: 2604 refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv] 2605 ParmDict: 2606 puts result F^2 in each ref[8] in refList 2607 ''' 2579 2608 global tsum 2580 2609 parmDict.update(dict(zip(varyList,values))) 2581 2610 Tdata,Xdata = GetAtomTX(RBdata,parmDict) 2582 2611 Mdata = parmDict['Mdata'] 2583 FF = np.zeros(len(Tdata))2584 2612 MDval = parmDict['0:MDval'] 2585 2613 MDaxis = parmDict['0:MDaxis'] … … 2587 2615 Srefs = np.zeros(len(refList)) 2588 2616 sumFcsq = 0. 2589 t0 = time.time()2590 2617 for refl in refList: 2591 for i,El in enumerate(Tdata): #NB: generator here is slower! 2592 FF[i] = refl[7][El] 2593 FF *= Mdata/len(refl[8]) #FF*occ 2594 phase = np.inner(refl[8],Xdata) #hx+ky+lz 2595 phase += refl[9][:,np.newaxis] #hx+ky+lz+phi 2596 cosp = np.cos(twopi*phase) 2597 fas = np.sum(FF*cosp) 2598 if ifInv: 2599 fbs = 0. 2600 else: 2601 sinp = np.sin(twopi*phase) 2602 fbs = np.sum(FF*sinp) 2603 refl[5] = (fas**2+fbs**2)*refl[3] #*calcMDcorr(MDval,MDaxis,Uniq,Gmat) 2618 t0 = time.time() 2619 refl[5] = mcsasfCalc(ifInv,Tdata,Mdata,Xdata,refl[3],refl[7],refl[8],refl[9]) 2620 # refl[5] *= calcMDcorr(MDval,MDaxis,Uniq,Gmat) 2621 tsum += (time.time()t0) 2604 2622 sumFcsq += refl[5] 2605 tsum += (time.time()t0)2606 2623 scale = (parmDict['sumFosq']/sumFcsq) 2607 2624 for iref,refl in enumerate(refList): … … 2654 2671 parmDict['atNo'] = atNo #total no. of atoms 2655 2672 parmDict['nObj'] = len(MCSAObjs) 2673 aTypes = list(aTypes) 2656 2674 Tdata,Xdata = GetAtomTX(RBdata,parmDict) 2657 2675 parmDict['Mdata'] = GetAtomM(Xdata,SGData) … … 2666 2684 SQ = 0.25/d**2 2667 2685 Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:] 2668 FFs = G2el.getFFvalues(FFtables,SQ )2686 FFs = G2el.getFFvalues(FFtables,SQ,True) 2669 2687 refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi]) 2670 2688 sumFosq += f*m … … 2692 2710 SQ = 0.25/d**2 2693 2711 Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:] 2694 FFs = G2el.getFFvalues(FFtables,SQ )2712 FFs = G2el.getFFvalues(FFtables,SQ,True) 2695 2713 refs.append([h,k,l,m,f*m,iref,0.,FFs,Uniq,phi]) 2696 2714 sumFosq += f*m … … 2724 2742 SQ = 0.25/d**2 2725 2743 Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:] 2726 FFs = G2el.getFFvalues(FFtables,SQ )2744 FFs = G2el.getFFvalues(FFtables,SQ,True) 2727 2745 refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi]) 2728 2746 sumFosq += f*m 2747 nRef = len(refs) 2729 2748 rcov = np.identity(len(refs)) 2749 print ' Minimum dspacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef) 2730 2750 parmDict['sumFosq'] = sumFosq 2731 2751 x0 = [parmDict[val] for val in varyList]
Note: See TracChangeset
for help on using the changeset viewer.