<commit_before>#!/opt/local/bin/python
u"""
Calculate the cohesive energy as a function of lattice constant,
by altering the lattice constant in pmd00000 file.

And if possible, calculate equilibrium lattice size and
bulk modulus, too.
"""

import sys,os,commands
import numpy as np
import optparse
from scipy.optimize import leastsq
import matplotlib.pyplot as plt


def read_pmd():
    f=open('0000/pmd00000','r')
    #...read 1st line and get current lattice size
    al= float(f.readline().split()[0])
    hmat= np.zeros((3,3))
    hmat[0]= [ float(x) for x in f.readline().split() ]
    hmat[1]= [ float(x) for x in f.readline().split() ]
    hmat[2]= [ float(x) for x in f.readline().split() ]
    f.readline()
    f.readline()
    f.readline()
    natm= int(f.readline().split()[0])
    f.close()
    return (al,hmat,natm)

def get_vol(al,hmat):
    a1= hmat[0:3,0] *al
    a2= hmat[0:3,1] *al
    a3= hmat[0:3,2] *al
    return np.dot(a1,np.cross(a2,a3))

def replace_1st_line(x):
    f=open('0000/pmd00000','r')
    ini= f.readlines()
    f.close()
    g=open('0000/pmd00000','w')
    for l in range(len(ini)):
        if l == 0:
            g.write(' {0:10.4f}\n'.format(x))
        else:
            g.write(ini[l])
    g.close()

def residuals(p,y,x):
    b,bp,v0,ev0= p
    err= y -( b*x/(bp*(bp-1.0)) *(bp*(1.0-v0/x) +(v0/x)**bp -1.0) +ev0 )
    return err

def peval(x,p):
    b,bp,v0,ev0= p
    return b*x/(bp*(bp-1.0)) *(bp*(1.0-v0/x) +(v0/x)**bp -1.0) +ev0

if __name__ == '__main__':
    
    usage= '%prog [options] <min> <max>'

    parser= optparse.OptionParser(usage=usage)
    parser.add_option("-n",dest="niter",type="int",default=10,
                      help="Number of points to be calculated.")
    parser.add_option("-p",action="store_true",
                      dest="plot",default=False,
                      help="Plot a graph on the screen.")
    parser.add_option("--pmdexec",dest="pmdexec",type="string",
                      default='~/src/nap/pmd/pmd',
                      help="path to the pmd executable.")
    (options,args)= parser.parse_args()

    if len(args) != 2:
        print ' [Error] number of arguments wrong !!!'
        print usage
        sys.exit()

    niter= options.niter
    shows_graph= options.plot
    pmdexec= options.pmdexec

    al_orig,hmat,natm= read_pmd()
    al_min = float(sys.argv[1])
    al_max = float(sys.argv[2])

    if al_orig < al_min or al_orig > al_max:
        print ' [Warning] min and max maybe wrong.'
        print '   hoping you know what you are doing.'
        print '   al_min, al_orig, al_max=',al_min, al_orig, al_max
        #sys.exit()

    logfile= open('log.Ecoh-vs-size','w')
    outfile1= open('out.Ecoh-vs-size','w')
    dl= (al_max -al_min)/niter
    for iter in range(niter+1):
        al= al_min +dl*iter
        replace_1st_line(al)
        os.system(pmdexec +' > out.pmd')<commit_msg>fatal bug in Ecoh-vs-size.py fixed.<commit_after>#!/opt/local/bin/python
u"""
Calculate the cohesive energy as a function of lattice constant,
by altering the lattice constant in pmd00000 file.

And if possible, calculate equilibrium lattice size and
bulk modulus, too.
"""

import sys,os,commands
import numpy as np
import optparse
from scipy.optimize import leastsq
import matplotlib.pyplot as plt


def read_pmd():
    f=open('0000/pmd00000','r')
    #...read 1st line and get current lattice size
    al= float(f.readline().split()[0])
    hmat= np.zeros((3,3))
    hmat[0]= [ float(x) for x in f.readline().split() ]
    hmat[1]= [ float(x) for x in f.readline().split() ]
    hmat[2]= [ float(x) for x in f.readline().split() ]
    f.readline()
    f.readline()
    f.readline()
    natm= int(f.readline().split()[0])
    f.close()
    return (al,hmat,natm)

def get_vol(al,hmat):
    a1= hmat[0:3,0] *al
    a2= hmat[0:3,1] *al
    a3= hmat[0:3,2] *al
    return np.dot(a1,np.cross(a2,a3))

def replace_1st_line(x):
    f=open('0000/pmd00000','r')
    ini= f.readlines()
    f.close()
    g=open('0000/pmd00000','w')
    for l in range(len(ini)):
        if l == 0:
            g.write(' {0:10.4f}\n'.format(x))
        else:
            g.write(ini[l])
    g.close()

def residuals(p,y,x):
    b,bp,v0,ev0= p
    err= y -( b*x/(bp*(bp-1.0)) *(bp*(1.0-v0/x) +(v0/x)**bp -1.0) +ev0 )
    return err

def peval(x,p):
    b,bp,v0,ev0= p
    return b*x/(bp*(bp-1.0)) *(bp*(1.0-v0/x) +(v0/x)**bp -1.0) +ev0

if __name__ == '__main__':
    
    usage= '%prog [options] <min> <max>'

    parser= optparse.OptionParser(usage=usage)
    parser.add_option("-n",dest="niter",type="int",default=10,
                      help="Number of points to be calculated.")
    parser.add_option("-p",action="store_true",
                      dest="plot",default=False,
                      help="Plot a graph on the screen.")
    parser.add_option("--pmdexec",dest="pmdexec",type="string",
                      default='~/src/nap/pmd/pmd',
                      help="path to the pmd executable.")
    (options,args)= parser.parse_args()

    if len(args) != 2:
        print ' [Error] number of arguments wrong !!!'
        print usage
        sys.exit()

    niter= options.niter
    shows_graph= options.plot
    pmdexec= options.pmdexec

    al_orig,hmat,natm= read_pmd()
    al_min = float(args[0])
    al_max = float(args[1])

    if al_orig < al_min or al_orig > al_max:
        print ' [Warning] min and max maybe wrong.'
        print '   hoping you know what you are doing.'
        print '   al_min, al_orig, al_max=',al_min, al_orig, al_max
        #sys.exit()

    logfile= open('log.Ecoh-vs-size','w')
    outfile1= open('out.Ecoh-vs-size','w')
    dl= (al_max -al_min)/niter
    for iter in range(niter+1):
        al= al_min +dl*iter
        replace_1st_line(al)
        os.system(pmdexec +' > out.pmd')
{'fraction_non_alphanumeric': 0.1692519820751465, 'fraction_numerical': 0.03136849362288866, 'mean_word_length': 2.615576323987539, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5490936', 'n_tokens_mistral': 2310, 'n_tokens_neox': 2203, 'n_words': 495}