dacase / nabc

Nonperiodic force field calculations, with a "C" language interface rather than the original interface based on the Nucleic Acid Builder (NAB) molecular manipulation language
Other
6 stars 1 forks source link

Some users may need a larger maximum system size in nmode.c #3

Open Nate-Avish opened 7 months ago

Nate-Avish commented 7 months ago

Hello,

I'm computing normal modes for a system with about 13,300 atoms. However, this system size overflows the INT_T "lwork" variable in nmode.c.

I've gotten around this issue by modifying a few files in the source code, but it isn't the prettiest fix, and other users with large systems may wish for an easier time than I had :)

Thank you for creating and maintaining this library! The newton() function and the portions of nmode() I've used have proven quite capable for my system.

Best wishes, Nate

dacase commented 7 months ago

On Thu, Feb 08, 2024, Nate-Avish wrote:

I'm computing normal modes for a system with about 13,300 atoms. However, this system size overflows the INT_T "lwork" variable in nmode.c.

I've gotten around this issue by modifying a few files in the source code, but it isn't the prettiest fix, and other users with large systems may wish for an easier time than I had :)

If you want to send along your changes, I'll try to put them in. When these codes were written, the notion that one could do normal modes on 13,000 atom systems seemed remote--computers did not generally have 20 Gb or more of main memory.

....dac

Nate-Avish commented 7 months ago

I can't send the exact files, since I've added a large number of print statements to some, which would make them more difficult to search. So I've listed the modifications below, including context. Also, please note I have next to no proficiency in C or Fortran, so my changes are not elegant, and I kept (commented out) the original lines as a sort of panic button :)

nabc/src/sff/nmode.c In ntrun == 1 section,

New line: long long int longwork = 1 + 6 * ncopy + 2 * ncopy * ncopy;

Context:

      if (eigp)
         jobz = 'V';
      else
         jobz = 'N';
      uplo = 'U';
      lwork = 1 + 6 * ncopy + 2 * ncopy * ncopy;                                               
      long long int longwork = 1 + 6 * ncopy + 2 * ncopy * ncopy;                                                                          
      liwork = 3 + 5 * ncopy;

New or commented out lines:

//work = vector(0, lwork);
work = longvector(0, longwork);

Context:

      liwork = 3 + 5 * ncopy;                                                                        
      //work = vector(0, lwork);                                              
      work = longvector(0, longwork);                                                                                                                                      
      iwork = ivector(0, liwork);

nabc/src/sff/memutil.c (and memutil.h, as appropriate)

New functions:

REAL_T  *longvector( size_t nl, long long int nh )
{
  REAL_T  *v;

  v = ( REAL_T * )malloc( ( nh - nl + 1 ) * sizeof( REAL_T ) );
  if( !v ) nrerror( "allocation failure in vector()" );
  return( v - nl );
}
void    free_longvector( REAL_T *v, size_t nl, long long int nh )
{
  free( v + nl );
}

nabc/src/lapack/dsyev.f

New line:

INTEGER(kind=SELECTED_INT_KIND(10)) :: LNGWRK

Context:

*     .. Local Scalars ..                                                                                 
      LOGICAL            LOWER, LQUERY, WANTZ
      INTEGER            IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
     $                   LLWORK, LWKOPT, NB

      INTEGER(kind=SELECTED_INT_KIND(10)) :: LNGWRK
      DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
     $                   SMLNUM
*     ..                                                                                                  
*     .. External Functions ..

New or commented out lines:

*      LLWORK = LWORK - INDWRK + 1                                                                        
      LNGWRK = LWORK - INDWRK + 1
*      CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),                                     
*     $             WORK( INDWRK ), LLWORK, IINFO )                                                       
      CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
     $             WORK( INDWRK ), LNGWRK, IINFO )

Context:

Search the lines I commented out.

nabc/src/lapack/dsytrd.f

New or commented out lines:

*      INTEGER            INFO, LDA, LWORK, N                                                             
      INTEGER            INFO, LDA, N

      INTEGER(kind=SELECTED_INT_KIND(10)) :: LWORK

Context: Search the line I commented out

Best wishes, Nate