Changeset 965


Ignore:
Timestamp:
Jun 24, 2013 12:20:48 PM (8 years ago)
Author:
vondreele
Message:

MCSA with fortran

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r962 r965  
    3232import GSASIIspc as G2spc
    3333import numpy.fft as fft
     34import pypowder as pyd
    3435
    3536sind = lambda x: np.sin(x*np.pi/180.)
     
    23172318        x0 = schedule.getstart_temp(best_state)
    23182319    else:
    2319         x0 = random.uniform(size=len(x0))*(upper-lower) + lower
     2320#        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
    23202321        best_state.x = None
    23212322        best_state.cost = numpy.Inf
     
    25062507            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
    25072508        return np.array(Mdata)
    2508 
     2509       
    25092510    def GetAtomTX(RBdata,parmDict):
    25102511        'Needs a doc string'
     
    25192520                parm = ':'+str(iatm)+key
    25202521                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]
    25222526        iatm = nfixAt
    25232527        for iObj in range(parmDict['nObj']):
     
    25262530                if parmDict[pfx+'Type'] == 'Vector':
    25272531                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
    2528                     aTypes = RBRes['rbTypes']
    25292532                    vecs = RBRes['rbVect']
    25302533                    mags = RBRes['VectMag']
     
    25342537                elif parmDict[pfx+'Type'] == 'Residue':
    25352538                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
    2536                     aTypes = RBRes['rbTypes']
    25372539                    Cart = np.array(RBRes['rbXYZ'])
    25382540                    for itor,seq in enumerate(RBRes['rbSeq']):
     
    25492551                    for j in range(3):
    25502552                        Xdata[j][iatm] = X[j]
    2551                     Tdata[iatm] = aTypes[i]
     2553                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
    25522554                    iatm += 1
    25532555            elif parmDict[pfx+'Type'] == 'Atom':
     
    25562558                    parm = pfx+key[1:]              #remove extra ':'
    25572559                    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]
    25592564                iatm += 1
    25602565            else:
     
    25632568   
    25642569    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       
    25712600       
    25722601    def mcsaCalc(values,refList,rcov,ifInv,RBdata,varyList,parmDict):
    2573 #        ''' Compute structure factors for all h,k,l for phase
    2574 #        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 refList
    2578 #        '''       
     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        '''       
    25792608        global tsum
    25802609        parmDict.update(dict(zip(varyList,values)))
    25812610        Tdata,Xdata = GetAtomTX(RBdata,parmDict)
    25822611        Mdata = parmDict['Mdata']
    2583         FF = np.zeros(len(Tdata))
    25842612        MDval = parmDict['0:MDval']
    25852613        MDaxis = parmDict['0:MDaxis']
     
    25872615        Srefs = np.zeros(len(refList))
    25882616        sumFcsq = 0.
    2589         t0 = time.time()
    25902617        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)
    26042622            sumFcsq += refl[5]
    2605         tsum += (time.time()-t0)
    26062623        scale = (parmDict['sumFosq']/sumFcsq)
    26072624        for iref,refl in enumerate(refList):
     
    26542671    parmDict['atNo'] = atNo                 #total no. of atoms
    26552672    parmDict['nObj'] = len(MCSAObjs)
     2673    aTypes = list(aTypes)
    26562674    Tdata,Xdata = GetAtomTX(RBdata,parmDict)
    26572675    parmDict['Mdata'] = GetAtomM(Xdata,SGData)
     
    26662684                SQ = 0.25/d**2
    26672685                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
    2668                 FFs = G2el.getFFvalues(FFtables,SQ)
     2686                FFs = G2el.getFFvalues(FFtables,SQ,True)
    26692687                refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi])
    26702688                sumFosq += f*m
     
    26922710                SQ = 0.25/d**2
    26932711                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
    2694                 FFs = G2el.getFFvalues(FFtables,SQ)
     2712                FFs = G2el.getFFvalues(FFtables,SQ,True)
    26952713                refs.append([h,k,l,m,f*m,iref,0.,FFs,Uniq,phi])
    26962714                sumFosq += f*m
     
    27242742                SQ = 0.25/d**2
    27252743                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
    2726                 FFs = G2el.getFFvalues(FFtables,SQ)
     2744                FFs = G2el.getFFvalues(FFtables,SQ,True)
    27272745                refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi])
    27282746                sumFosq += f*m
     2747        nRef = len(refs)
    27292748        rcov = np.identity(len(refs))
     2749    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
    27302750    parmDict['sumFosq'] = sumFosq
    27312751    x0 = [parmDict[val] for val in varyList]
Note: See TracChangeset for help on using the changeset viewer.