Re: Bug in NAMD

From: Pawel Kedzierski (pawel.kedzierski_at_pwr.edu.pl)
Date: Fri Jun 21 2019 - 07:01:03 CDT

W dniu 20.06.2019 o 16:46, Marcelo C. R. Melo pisze:
> 2) How are you checking the coordinates sent to the QM software? Are
> you saving a DCD at every step?

I have confirmed that the coordinates saved by NAMD in the QM-only DCD
file are repeated every few steps.

Python code used to check this is pasted below.

> 3) I noticed you are using the NAMD build from Nov/2018. (NAMD 2.13
> for Linux-x86_64-multicore, Built Fri Nov 9 14:39:16 CST 2018 by jim
> on ganymede.ks.uiuc.edu <http://ganymede.ks.uiuc.edu/>)
> Do you have access to a new one? Could you download the pre-compiled
> one from the NAMD website for tests?

Downloaded and tested this binary build:

NAMD Git-2019-06-21 for Linux-x86_64-multicore
Based on Charm++/Converse 60800 for multicore-linux64-iccstatic
Built Fri Jun 21 02:03:40 CDT 2019 by jim on belfast.ks.uiuc.edu

Same problem.

Best regards,

Paweł Kędzierski

#!/usr/bin/env python3
# Load coordinates from DCD trajectory and compare them. Requires catdcd,
# Python >=3.5 and numpy Python library
# Arguments: file.psf file.dcd
import numpy as np, subprocess as sp, sys

inpsf, indcd = sys.argv[1:3]
sp.run(['catdcd', '-o', 'tmp.xyz', '-otype', 'xyz', '-s', inpsf,
         '-stype', 'psf', indcd], check=True, stdout=sp.PIPE)
prec = None
frames = []
with open('tmp.xyz', 'r') as xyz:
   for line in xyz:
     words = line.split()
     if len(words)==1: # the line with atom count
        natoms = int(words[0])
        next(xyz) # skip comment
        coords = np.zeros((natoms,3), dtype=float) # array for coordinates
        for i in range(natoms):
            words = next(xyz).split()
            if not prec: # recognize precision of coords in xyz file
               prec = len(words[1].split('.')[1])
            coords[i] = words[1:4] # read x,y,z line by line
        frames.append(coords) # remember this frame of coordinates
N = len(frames)
print("\nSuccessfully read %d frames with precision to the %d'th decimal point."
       % (N, prec))

# Now, compare the frames
eps = eval('1e-%d' % prec)
found = False
for i in range(N-1):
    anynew = False
    for j in range(i+1, N):
       if np.all(np.abs(frames[i] - frames[j]) < eps):
          if not found:
             print('Pairs of identical frames (0-based indices):')
             found = True
          anynew = True
          print('%d-%d' % (i,j), end=' ')
    if anynew:
       print() # flush line
if not found:
    print('None of the frames were identical.')


This archive was generated by hypermail 2.1.6 : Tue Sep 17 2019 - 23:20:52 CDT