MaginnGroup / PyLAT

GNU General Public License v3.0
101 stars 54 forks source link

GKC calculation on MD trajectory with out partial charge definition #12

Closed pythonpanda2 closed 3 years ago

pythonpanda2 commented 3 years ago

Hi,

I am running a ML forcefield to generate MD trajctory for computing GKC. There is no partial charge definition in input to LAMMPS simulation. I have written the trajectory in metal units in the following format

ITEM: ATOMS id type element x y z vx vy vz fx fy fz

I read in a closed post that 'id' can be changed to 'mol' in the PyLat line. I can also do the unit conversion in a post processing step.

I have tried to modify the lammps input data file by adding a column for partial charges. But this did not seem to work out. I am interested to know how can a charge array be passed to GKC calculation. Please let me know if you have any advice on this.

pythonpanda2 commented 3 years ago

This is the error that i get when i try to get the charges array can be seen below

  (molcharges, atomcharges,n) = gc.getmolcharges(datfilename,n)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-14-04db2081846f> in <module>
      2 output['Conductivity']['units'] = 'S/m'
      3 n = gc.findnumatoms(datfilename)
----> 4 (molcharges, atomcharges,n) = gc.getmolcharges(datfilename,n)

~/PyLAT_MOD/src/getAtomCharges.py in getmolcharges(self, datfilename, n)
     72         molcharges = np.zeros(nummol)
     73         for atom in range(0,n):
---> 74             molcharges[int(mol[int(atom)])-1] += atomcharges[int(atom)]
     75
     76         datfile.close()

IndexError: index -1 is out of bounds for axis 0 with size 0

The datafile looks like below

   1024  atoms
       2  atom types

  0.00000000      29.39400839  xlo xhi
  0.00000000      29.39400839  ylo yhi
  0.00000000      29.39400839  zlo zhi

Masses

       1   6.94000000    # Li
       2   35.45000000    # Cl

Atoms # charge

     1    1   1.000000        1.56564969       3.26370015       2.54800308
     2    2  -1.000000        3.71432217       4.03406815       3.01543742
     3    2  -1.000000        1.63489215       1.26760213       3.83945239
     4    2  -1.000000        1.35354083       2.68599049       0.08215336
mike5603 commented 3 years ago

How many blank lines are between the Atoms line and when the data begins? PyLAT assumes that there is exactly one blank line between them

mike5603 commented 3 years ago

The formatting on github confused me and there is only one blank line. Is this system atomic or molecular?

pythonpanda2 commented 3 years ago

data.txt

Thanks a lot for your feedback. I have attached the input file.

pythonpanda2 commented 3 years ago

The system is an LiCl melt

mike5603 commented 3 years ago

Since you have an atomistic system, add the following lines below line 66 (mol[int(line[0])-1] = int(line[1])) in getAtomCharges.py:

        elif len(line) == 9 or len(line) == 6:
            atomcharges[int(line[0])-1] = float(line[2])
            mol[int(line[0])-1] = int(line[0])

Let me know if this works

pythonpanda2 commented 3 years ago

Many thanks for suggesting the changes. I am running a test now. I will report back as soon as the test is finished.

pythonpanda2 commented 3 years ago

The COM runs finished. This time got a new error :(

COM velocity calculation 100.00% complete
begining GK conductivity correlation

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-13-11280baa795a> in <module>
----> 1 output = cc.calcConductivity(molcharges, trjfilename, logfilename, datfilename, temp, output, moltype, moltypel,verb,GKC_skip,GKC_Tolerance,GKC_J_Output)

~/PyLAT_MOD/src/calcCond.py in calcConductivity(self, molcharges, trjfilename, logfilename, datfilename, T, output, moltype, moltypel, ver, firstpoint, tol, Jout)
     66         if ver >= 2:
     67             print('begining GK conductivity correlation')
---> 68         J = self.clacJ(jx, jy, jz, dt, tsjump, firstpoint,ver)
     69         if Jout:
     70             self.writeJ(J,tsjump,dt)

~/PyLAT_MOD/src/calcCond.py in clacJ(self, jx, jy, jz, dt, tsjump, firstpoint, ver)
    193         for l in range(0,len(jx)):
    194             for m in range(0, len(jx)):
--> 195                 Jtest = self.correlate(jx[l],jx[m])
    196                 J[0] += Jtest
    197                 J[m+1] += Jtest

~/PyLAT_MOD/src/calcCond.py in correlate(self, a, b)
    350         bl = np.concatenate((b,np.zeros(len(b))),axis=0)
    351         c= np.fft.ifft(np.fft.fft(al)*np.conjugate(np.fft.fft(bl))).real
--> 352         d = c[:len(c)/2]
    353         d/=(np.arange(len(d))+1)[::-1]
    354         return d

TypeError: slice indices must be integers or None or have an __index__ method
mike5603 commented 3 years ago

That looks to be a relic of this code being written for python 2. I pushed a fix to master, but if you want to keep the changes, at the end of the calcCond.py code, replace d = c[:len(c)/2] with d = c[:int(len(c)/2)]

pythonpanda2 commented 3 years ago

Thanks a lot. The last fix worked out well. I have managed to finish one full batch of test.

jnitlion commented 3 years ago

Hi! I am having the same problem while trying to perform GKC calculations. The error says:

Traceback (most recent call last): File "/path/PyLAT.py", line 246, in (molcharges, atomcharges,n) = gc.getmolcharges(datfilename,n) File "/path/PyLAT-master/src/getAtomCharges.py", line 77, in getmolcharges molcharges[int(mol[int(atom)])-1] += atomcharges[int(atom)] IndexError: index -1 is out of bounds for axis 0 with size 0

which I think is same as the error raised in this issue. I have a coarse-grained system (spheres) with charges and dipole moment. Attached are my trajectory and input files for your reference. Also, is there a specific values or info needed inside the log file?

I hope you can help me. Thanks!

trajectory.txt data.txt