/aquad/src/AquaFuncts.c

Large files files are truncated, but you can click here to view the full file

/*
**  File: AquaFuncts.c
**
**  This file is part of the AQUA program suite, which is being
**  developed in the NMR Spectroscopy Department, Bijvoet Center for 
**  Biomolecular Research, Utrecht. 
**
**
**  AUTHOR(S) : Ton Rullmann
**  VERSION   : 1.5.0
**  DATE      : Dec 18, 1997
**
**  This file contains miscellaneous functions dealing with 
**  among other things Aqua data types.
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include "Compiler.h"
#include "General.h"

#include "AquaMacros.h"
#include "AquaTypes.h"
#include "AquaData.h"

/****************  function prototypes  ****************/

#include "AquaFuncts_cv.h"
#include "AquaFuncts_io.h"
#include "AquaStrucset.h"
#include "AquaFiles.h"
#include "MenuFuncts.h"
#include "Functs.h"

#include "AquaFuncts.h"

static Atom *GuessAtom( Atom *at );
static int CmpResidOfAtomptr( const void *a, const void *b );
static Mol *ExpandAtoms( Atom *atoms, const int at_count, Mol *molptr );
static Mol *MakeDefsMol( Atom *new_atoms, const int new_count, const Boolean makeres );

/****************  globals  ****************/

extern Boolean  iaf;                   /* import from Functs */
extern int      terseness;             /* import from Functs */
extern Card     msg_string;            /* import from Functs */

char            *project_name;         /* export to main program and libs */

#define max_MATCH 10
static Atomptr  MatchList[max_MATCH];
static int      find_atom_err;         /* flag FindAtom2 -> WriteNotFoundAtom */

/* BEGIN_EXTH */

/****************  Check_General_Options  ****************/
/*
**  Check the validity of the values of some general options against
**  the allowed range
*/

Boolean Check_General_Options ( void )
{
    
    Boolean check = true;    

    /* Terseness value [0-9] */
     if ( ( GetOptionIVal( opt_terseness ) < 0 ) || 
          ( GetOptionIVal( opt_terseness ) > 9 ) ) {
	check = false;
	printf("Option opt_terseness (%d) is outside it's bounds\n",GetOptionIVal(opt_terseness));
    }
   
    return ( check );
}

/****************  ProjectInit  ****************/
/*
**  Get the project name from the environment variable AQUAPROJECT.
*/

void ProjectInit( void )
{
    project_name = GetVar( "AQUAPROJECT", false );
    if ( !project_name )
    {
        puts( " * Define a project name first!" );
        exit( 1 );
    }
    if ( terseness < 6 ) {    
    printf( "Project name: %s\n", project_name );
    }
    return;
}

/****************  AquaDataFile  ****************/
/*
**  Open Aqua data file <file> in the AQUADATADIR directory.
*/

FILE *AquaDataFile( char *file )
{
    FILE    *fp;
    Filnam  filnam;
    char    *dir;
/* the total length of the string should be checked against size_FILNAM */

    dir = GetVar( "AQUADATADIR", true );
    strcpy( filnam, dir );
    strcat( filnam, DIRSEP );
    strcat( filnam, file );

    if ( terseness < 6 ) {
    printf( "Now opening: %s\n", filnam );
    }
    
    if ( NOTOPENr( fp, filnam ) )
    {
	printf( " * Aqua data-file <%s> could not be opened\n", filnam );
	if ( !iaf )
	    exit( 1 );
	return( NULL );
    }

    return( fp );
}

/****************  AquaFileName  ****************/
/*
**  Generate an Aqua file name.
**  The string is stored in <filnam> (which should contain at least
**  size_FILNAM characters), and a pointer to it is returned.
*/

char *AquaFileName( const Qtyp type, Typnam ident, char *filnam )
{
/* the total length of the string should be checked against size_FILNAM */

    strcpy( filnam, project_name );
    strcat( filnam, "_" );
    strcat( filnam, QtypSel( type ) );
    strcat( filnam, "_" );
    strcat( filnam, ident );
    strcat( filnam, FiltypExt( QtypFil( type ) ) );
    
    return( filnam );
}

/****************  OpenAquaFile  ****************/
/*
**  Generate the name of an Aqua file of Qtyp <type> by inputting the "ident" 
**  (if <reuse> is false) or using the last input ident (if <reuse> is true) 
**  and calling AquaFileName, 
**  and open the file in mode <mode>.
**  The return value is the file pointer, or null.
**  <filnam> is the address of a string that receives the file name; it should
**  contain at least size_FILNAM characters. 
*/

FILE *OpenAquaFile( const Qtyp type, char *mode, char *filnam, const Boolean reuse )
{
    FILE           *fp;
    Menu           *oldmenu = NULL;
    static Typnam  ident;
    char           *filnam_temp;

    if ( !reuse )
    {
        oldmenu = SetOneMenu( "Give file identification string:" ); 
        MReadWord( ident );
    }
    
    filnam_temp = AquaFileName( type, ident, filnam );
    
    if ( terseness < 6 ) {
        if ( *mode == 'a' ) {
            printf( "Now (re)opening: %s\n", filnam_temp);
        } else {
            printf( "Now opening: %s\n", filnam_temp );
        }
    }
        
    if ( ( fp = fopen( filnam, mode ) ) == NULL )
    {
        printf( " * File <%s> could not be opened in mode <%s>\n", filnam, mode );
        if ( !iaf )
            exit( 1 );
        return( NULL );
    }
    if ( !reuse )
        ResetMenu( oldmenu );
    return( fp );
}

/****************  ReadAquaIDLines  ****************/
/*
**  Read the standard Aqua identification lines.
**  Returns true if successful, otherwise false.
*/

Boolean ReadAquaIDLines( FILE *fp )
{
/***
    Card  version;
    if ( fscanf( fp, "$ Aqua Version %s", version ) < 1 )
    {
	puts( " * This is not an Aqua file" );
	return( false );
    }
    if ( NOTEQ( version, aqua_version ) )
    {
	puts( " * This file was created by an old Aqua version" );
	return( false );
    }
***/

    float  version;
    if ( fscanf( fp, "$ Aqua Version %f", &version ) < 1 )
    {
	puts( " * This is not an Aqua file" );
	return( false );
    }
    if ( version < aqua_versnum_comp )
    {
	puts( " * This file was created by an old Aqua version" );
	return( false );
    }
    else if ( version > aqua_versnum )
    {
	puts( " * This file was created by a future Aqua version" );
	exit( 1 );
    }
    FlushLine( fp );
    FlushLine( fp );
    return( true );
}

/****************  WriteAquaIDLines  ****************/
/*
**  Write the standard Aqua identification lines
*/

void WriteAquaIDLines( FILE *fp, const Qtyp qtype, Filnam filnam )
{
    time_t   tim;

    tim = time( NULL );
    fprintf( fp, "$ Aqua Version %s DATA: %s\n", aqua_version, 
             QtypTxt( qtype ) );
    fprintf( fp, "$ FILE: %s CREATED: %s", filnam, ctime( &tim ) );
}

/****************  ProcessAquaLine  ****************/
/*
**  Read until a non-blank line and process as standard Aqua header line. 
**  Format assumed:
**  $ DESCRIPTOR: DATA
**  The line is input stored at <line>.
**  DESCRIPTOR and DATA are stored at <descr> and <inp>.
**  These 3 first arguments should be pointers to non-overlapping 
**  character arrays of at least <size> bytes.
**  The function returns NULL if an error occurs, otherwise <line>.
*/

char *ProcessAquaLine( char *line, char *descr, char *inp, 
                       const int size, FILE *fp )
{
    char  *colptr;

    do
	if ( !ReadLine( line, size, fp ) )
	{
	    puts( " * Error while reading" );
	    return( NULL );
	}
    while
	( IsEmptyStr( line ) );

    colptr = strchr( line, ':' );
    if ( *line != '$' || !colptr )
    {
	puts( " * Unexpected line encountered:" );
	puts( line );
	return( NULL );
    }
    strcpy( descr, FirstNonBlank( line+1 ) );
    strcpy( inp, FirstNonBlank( colptr+1 ) );
    if ( OptionTrue( opt_prtaq ) )
	printf( "    read: %s\n", line );
    return( line );
}

/****************  QtypSel  ****************/
/*
**  Return file selector string belonging to aqua type <type>
*/

char *QtypSel( const Qtyp type )
{
    unsigned int i;
    
    for ( i=0; i<size_aquarec; i++)
        if ( type == aquarec[i].type )
            return( aquarec[i].file_selector );
            
    return( aquarec[0].file_selector );
}

/****************  QtypTxt  ****************/
/*
**  Return type text belonging to aqua type <type>
*/

char *QtypTxt( const Qtyp type )
{
    unsigned int i;
    
    for ( i=0; i<size_aquarec; i++)
        if ( type == aquarec[i].type )
            return( aquarec[i].type_text );
            
    return( aquarec[0].type_text );
}

/****************  QtypFil  ****************/
/*
**  Return file type belonging to aqua type <type>
*/

Filtyp QtypFil( const Qtyp type )
{
    unsigned int i;

    for ( i=0; i<size_aquarec; i++)
        if ( type == aquarec[i].type )
            return( aquarec[i].file_type );
            
    return( aquarec[0].file_type );
}

/****************  FiltypNam  ****************/
/*
**  Return name belonging to file type <type>
*/

char *FiltypNam( const Filtyp type )
{
    unsigned int i;
    
    for ( i=0; i<size_filerec; i++)
        if ( type == filerec[i].type )
            return( filerec[i].type_name );
            
    return( filerec[0].type_name );
}

/****************  FiltypLIBNam  ****************/
/*
**  Return library name belonging to file type <type>
*/

char *FiltypLIBNam( const Filtyp type )
{
    switch ( type )
    {
	case pdb_type:
	    return( "pdb" );
	case pdbn_type:
	    return( "pdbn" );
	case pdbx_type:
	    return( "pdbx" );
	case pdbm_type:
	    return( "pdbm" );
	case biosym_type:
	    return( "biosym" );
	default:
	    return( "unknown" );
    }
}

/****************  FiltypExt  ****************/
/*
**  Return file name extension belonging to file type <type>
*/

char *FiltypExt( const Filtyp type )
{
    unsigned int i;
    
    for ( i=0; i<size_filerec; i++)
        if ( type == filerec[i].type )
            return( filerec[i].extension );
            
    return( filerec[0].extension );
}

/****************  ContactRangeCount  ****************/
/*
**  Return number of ContactRange types defined.
*/

int ContactRangeCount( void )
{
    return( size_contactrec );
}

/****************  ContactRangeNam  ****************/
/*
**  Return name belonging to ContactRange type <type>.
*/

char *ContactRangeNam( const ContactRange type )
{
    unsigned int i;
    
    for ( i=0; i<size_contactrec; i++)
        if ( type == contactrec[i].type )
            return( contactrec[i].type_name );
            
    return( contactrec[0].type_name );
}

/****************  ContactRangeStarNam  ****************/
/*
**  Return NMR-STAR name belonging to ContactRange type <type>.
*/

char *ContactRangeStarNam( const ContactRange type )
{
    unsigned int i;
    
    for ( i=0; i<size_contactrec; i++)
        if ( type == contactrec[i].type )
            return( contactrec[i].star_type_name );
            
    return( contactrec[0].star_type_name );
}

/****************  ContactRangeIndex  ****************/
/*
**  Return index belonging to ContactRange type <type>.
*/

int ContactRangeIndex( const ContactRange type )
{
    unsigned int i;
    
    for ( i=0; i<size_contactrec; i++)
        if ( type == contactrec[i].type )
            return( i );
            
    return( -1 );
}

/****************  MapAtomsInSet  ****************/
/*
**  For a strucset map the 'external' atom names to the 
**  'internal' (Aqua) ones, or vice versa (determined by <how>).
**  If needed, struc_mol's are created and a loop over them is performed.
**  Cf. MapAtomsInDRestraints.
*/

void MapAtomsInSet( Strucset *set, const int how, char *lib_nam )
{
    char       *fct_name = "MapAtomsInSet";
    Struc      *strucptr;
    Mol        *molptr = set->molecule;
    Filnam     predef_nam;

    sprintf( msg_string, "preparing to map atom names in strucset <%s>",
             set->name );
    if ( terseness < 6 ) {
        STATUSm;
    }

    if ( set->mapped )
    {
        if ( terseness < 6 ) {
	    STATUS( "atom names already mapped" );
        }
	return;
    }

/* predefine the library name for use with the "whatif" option */

    strcpy( predef_nam, project_name );
    strcat( predef_nam, FiltypExt( atlib_type ) );

/* DO THIS HERE ? >>>>>>>>>>>>>>>>>>>>>>>> */

    if ( EQn( lib_nam, "whatif", 6 ) && set->struc_count > 1 )
	set->onemap = false;

/* use the unique molecule */

    if ( set->onemap )
	MapAtomsInMol( molptr, how, lib_nam, predef_nam );

/* or use the model-dependent molecules */
/* THE MAP IN set->molptr LOSES ITS MEANING */
    else
    {
	strucptr = set->fst_struc;
	while ( strucptr )
	{
	    if ( !strucptr->struc_mol )       /* create if necessary */
	    {
		strucptr->struc_mol = MolCreate( set->atom_count );
		MolCopy( strucptr->struc_mol, molptr );
		AddAtomSToMol( strucptr->struc_mol, molptr->atoms, molptr->atom_count,
		               set->atom_count );
	    }
	    if ( set->struc_count != 1 )      /* name as given by whnamchk */
		sprintf( predef_nam, "%s_%03i%s", project_name, 
		         strucptr->struc_num, FiltypExt( atlib_type ) );
	    MapAtomsInMol( strucptr->struc_mol, how, lib_nam, predef_nam );
	    strucptr = strucptr->next_struc;
	}
    }

    set->mapped = true;
}

/****************  MapAtomsInMol  ****************/
/*
**  Map the 'external' atom names to the 'internal' (Aqua) ones, 
**  or vice versa (determined by <how>).
**  Cf. MapAtomsIn1DR.
**  The library name <lib_nam> is derived from the coordinate file type,  
**  if it is "". Names are just copied if <lib_nam> is "copy" or "aqua". 
**
**  Name guessing - as implemented in MapAtoms->FindAtom - is not performed.
*/

void MapAtomsInMol( Mol *molptr, const int how, char *lib_nam_orig, Filnam predef_nam )
{
    char       *fct_name = "MapAtomsInMol";
    Mol        *atmap_mol = NULL, *atmap_mol2 = NULL;
    FctFindAt  Find = NULL;
    Strucset   tmpSet;
    Card       lib_nam;
    char       *lib_nam1 = lib_nam, *lib_nam2;

/* make a list of atom name pairs from a library file */
/* the library name is derived from the coordinate file type, if it is "" */ 
/* if the file type is "aqua", names will just be copied */

/* was a combination of libraries specified? if so, split the name on the + character */

    strcpy( lib_nam, lib_nam_orig );
    if ( ( lib_nam2 = strchr( lib_nam, '+' ) ))
    {
	*lib_nam2 = '\0';      /* this sets lib_nam1 to a correct string */
	lib_nam2++;
    }
    
/*
    printf( "first  lib: <%s>\n", lib_nam1 );
    if ( lib_nam2 )
	printf( "second lib: <%s>\n", lib_nam2 );
    else
	puts( "no second lib" );
*/

    if ( NOTEQ( lib_nam1, "copy" ) && NOTEQ( lib_nam1, "aqua" ) )
    {
	if ( EQ( lib_nam1, "" ) )
	    lib_nam1 = FiltypLIBNam( molptr->inp_filtyp );
	atmap_mol = ProcessAtomLIB( molptr, lib_nam1, predef_nam );
    }

/* second library - handling same as in MapAtomsIn1DR */

    if ( lib_nam2 && atmap_mol )
    {
	if ( IsEmptyStr( lib_nam2 ) )
	    lib_nam2 = FiltypLIBNam( molptr->inp_filtyp );
/*	    lib_nam2 = GetLIBName( item_NONE );              */
	if ( ( atmap_mol2 = ProcessAtomLIB( molptr, lib_nam2, predef_nam ) ) )
	{
	    StrucsetInit( &tmpSet, "allDefs", false );
	    StrucsetCreate( &tmpSet, atmap_mol->atom_count + atmap_mol2->atom_count, 0 );
	    AddAtomSToSet( &tmpSet, atmap_mol->atoms, atmap_mol->atom_count );
	    AddAtomSToSet( &tmpSet, atmap_mol2->atoms, atmap_mol2->atom_count );
	    DeleteMol( atmap_mol2 );
	    if ( OptionTrue( opt_find2 ) )
		MakeResidues( &tmpSet, false );     /* cf. MakeDefsMol */
	    atmap_mol = tmpSet.molecule;
	    sprintf( msg_string, "%i atom name definitions active for present molecule", 
		 atmap_mol->atom_count );
            if ( terseness < 6 ) {
	        STATUSm;
            }
	    if ( OptionTrue( opt_prtlib ) )
		WriteFAtoms( stdout, atmap_mol->atoms, atmap_mol->atom_count );
	}
    }

/* print what is going to happen */

    if ( !atmap_mol )
    {
        if ( terseness < 6 ) {
	    if ( how == int_to_ext )
	        STATUS( "copying internal to external atom names" );
	    else
	        STATUS( "copying external to internal atom names" );
        }
    }
    else
    {
	if ( how == int_to_ext )
	{
	    if ( terseness < 6 ) {
	        STATUS( "mapping internal to external atom names" );
            }
	    if ( OptionTrue( opt_find2 ) )
	    {
		if ( terseness < 6 ) {
		    STATUS( "using routine FindAtom2" );
                }
		Find = FindAtom2;
	    }
	    else
	    {
		if ( terseness < 6 ) {
		    STATUS( "using routine FindAtom3" );
                }
		Find = FindAtom3;
	    }
	}
	else
	{
	    if ( terseness < 6 ) {
	        STATUS( "mapping external to internal atom names" );
            }
	    if ( OptionTrue( opt_find2 ) )
	    {
		if ( terseness < 6 ) {
		    STATUS( "using routine FindAtom2" );
                }
		Find = FindAtom2;
	    }
	    else
	    {
		if ( terseness < 6 ) {
		    STATUS( "using routine FindAtom3" );
                }
		Find = FindAtom3;
	    }
	}
    }

/* now do the mapping or copying */

    MapAtoms( molptr->atoms, molptr->atom_count, atmap_mol, how, false, Find );

/* free the list again */

    if ( atmap_mol )
	DeleteMol( atmap_mol );
}

/****************  MapAtoms  ****************/
/*
**  In <atoms> map the 'external' atom names to the 'internal' (Aqua) 
**  ones, or v.v., according to the names of the atoms appearing in <atmap_mol>.
*/

void MapAtoms( Atom *atoms, int atom_count, Mol *atmap_mol, 
               const int how, const Boolean guess, FctFindAt FindAtom )
{
    Atom     *atmap;
    Boolean  err = false;
    int      idx, find_which;

    if ( !atmap_mol )
    {

/* if no list was made, just copy */

	if ( how == int_to_ext )
	{
	    while ( atom_count-- )
	    {
		memcpy( atoms->atom_ext, atoms->atom_int, size_ID );
		atoms++;
	    }
	}
	else
	{
	    while ( atom_count-- )
	    {
		memcpy( atoms->atom_int, atoms->atom_ext, size_ID );
		atoms++;
	    }
	}
	return;
    }


/* if a list was made, use it to map */

    find_which = ( how == int_to_ext ? use_int : use_ext );
    while ( atom_count-- )
    {
	if ( !( atmap = FindAtom( &idx, atoms, atmap_mol, find_which, guess ) ) )
	{
	    WriteNotFoundAtom( stdout, atoms, find_which, find_which );
	    puts( "copying atom name instead" );
	    if ( how == int_to_ext )
		memcpy( atoms->atom_ext, atoms->atom_int, size_ID );
	    else
		memcpy( atoms->atom_int, atoms->atom_ext, size_ID );
	    err = true;
	}
	else
	{
	    if ( how == int_to_ext )
		memcpy( atoms->atom_ext, atmap->atom_ext, size_ID );
	    else
		memcpy( atoms->atom_int, atmap->atom_int, size_ID );
	}
	atoms++;
    }
    
    if ( err )
	puts( " * One or more atoms could not be found in the AtomLIB" );

}

/****************  MapPtrAtoms  ****************/
/*
**  In <atomptrs> map the 'external' atom names to the 'internal' (Aqua) 
**  ones, or v.v., according to the names of the atoms appearing in <atmap_mol>.
*/

void MapPtrAtoms( Atomptr *atomptrs, int atom_count, Mol *atmap_mol, 
                  const int how, const Boolean guess, FctFindAt FindAtom )
{
    Atom     *atmap;
    Boolean  err = false;
    int      idx, find_which;

    if ( !atmap_mol )
    {

/* if no list was made, just copy */

	if ( how == int_to_ext )
	{
	    while ( atom_count-- )
	    {
		memcpy( (*atomptrs)->atom_ext, (*atomptrs)->atom_int, size_ID );
		atomptrs++;
	    }
	}
	else
	{
	    while ( atom_count-- )
	    {
		memcpy( (*atomptrs)->atom_int, (*atomptrs)->atom_ext, size_ID );
		atomptrs++;
	    }
	}
	return;
    }


/* if a list was made, use it to map */

    find_which = ( how == int_to_ext ? use_int : use_ext );
    while ( atom_count-- )
    {
	if ( !( atmap = FindAtom( &idx, *atomptrs, atmap_mol, find_which, guess ) ) )
	{
	    WriteNotFoundAtom( stdout, *atomptrs, find_which, find_which );
	    puts( "copying atom name instead" );
	    if ( how == int_to_ext )
		memcpy( (*atomptrs)->atom_ext, (*atomptrs)->atom_int, size_ID );
	    else
		memcpy( (*atomptrs)->atom_int, (*atomptrs)->atom_ext, size_ID );
	    err = true;
	}
	else
	{
	    if ( how == int_to_ext )
		memcpy( (*atomptrs)->atom_ext, atmap->atom_ext, size_ID );
	    else
		memcpy( (*atomptrs)->atom_int, atmap->atom_int, size_ID );
	}
	atomptrs++;
    }
    
    if ( err )
	puts( " * One or more atoms could not be found in the AtomLIB" );

}

/****************  MakeResidues  ****************/
/*
**  
*/

void MakeResidues( Strucset *set, Boolean sort )
{
    char      *fct_name = "MakeResidues";
    Mol       *molptr;
    Atom      *atoms;
    Atomptr   *atomptrs, *list;
    Residue   *residues;
    int       *atomidxs, *list2;
    int       i, j, idx, res_count, atom_count;
    char      *res_id;

    molptr = set->molecule;
    atom_count = molptr->atom_count;

/* make a list of pointers to the atoms */

    list = CREATE_NEW( Atomptr, atom_count );
    if ( !list )
	NOMEM( "atomptrs" );
    atomptrs = list;
    atoms = molptr->atoms;
    for ( i=0; i < atom_count; i++ )
	*atomptrs++ = atoms++;

/* allocate a list of atom indices */

    list2 = CREATE_NEW( int, atom_count );
    if ( !list2 )
	NOMEM( "atomidxs" );

/* make a unique residue ID name */

    atoms = molptr->atoms;
    for ( i=0; i < atom_count; i++ )
	MakeUniqResid( atoms++ );

/* order the list on unique residue IDs */

    if ( sort )
	qsort( list, (size_t)atom_count, sizeof( Atomptr ), CmpResidOfAtomptr );

/* now count the number of residues and set the new, sequential residue number */
/* fill the atom index array */

    atomidxs = list2;
    atomptrs = list;
    atoms = *atomptrs;

    res_id = atoms->res_id;
    res_count = 1;

    for ( i=0; i < atom_count; i++ )
    {
        if ( NOTEQ( atoms->res_id, res_id ) )
        {
            res_id = atoms->res_id;
            res_count++;
        }
	atoms->res_num_new = res_count;   /* sequential number, starting at 1 */
	for ( j=0; j < atom_count; j++ )
	{
	    if ( atoms == &molptr->atoms[j] )
	    {
		*atomidxs++ = j;
		break;
	    }
	}
        atoms = *(++atomptrs);
    }

/* make the residue list */

    sprintf( msg_string, "allocating %i residues in molecule of strucset <%s>",
             res_count, set->name );
    if ( terseness < 6 ) {
        STATUSm;
    }
    residues = CREATE_NEW( Residue, res_count );
    
    molptr->residues = residues;
    molptr->res_count = res_count;

    atomidxs = list2;
    atomptrs = list;
    atoms = *atomptrs;

    res_id = atoms->res_id;
    residues->res_num_new  = atoms->res_num_new;
    residues->res_num_orig = atoms->res_num_orig;
    memcpy( residues->chain,   atoms->chain,   size_ID );
    memcpy( residues->res_nam, atoms->res_nam, size_ID );
    strcpy( residues->res_id,  atoms->res_id );
    residues->res_insert = atoms->res_insert;
    residues->atomptrs = atomptrs;
    residues->atomidxs = atomidxs;

    for ( i=0; i < atom_count; i++ )
    {
        if ( NOTEQ( atoms->res_id, res_id ) )
        {
            residues++;

            res_id = atoms->res_id;
            residues->res_num_new  = atoms->res_num_new; 
            residues->res_num_orig = atoms->res_num_orig;
            memcpy( residues->chain,   atoms->chain,   size_ID );
            memcpy( residues->res_nam, atoms->res_nam, size_ID );
            strcpy( residues->res_id,  atoms->res_id );
            residues->res_insert = atoms->res_insert;
            residues->atomptrs = atomptrs;
            residues->atomidxs = atomidxs;
        }
        residues->atom_count++;
        atoms = *(++atomptrs);
	++atomidxs;
    }

/* check */

    atoms = molptr->atoms;
    atomidxs = list2;
    atomptrs = list;
    for ( i=0; i < atom_count; i++ )
    {
	idx = *atomidxs++;
	if ( ( atoms + idx ) != *atomptrs++ )
	{
	    ERROR( " * Atom index incorrect - code error" );
	    exit( 1 );
	}	    
    }

/* print */

    if ( OptionTrue( opt_prtres ) )
        WriteResidues( stdout, molptr->residues, molptr->res_count );
}

/****************  MakeUniqResid  ****************/
/*
**  
*/

void MakeUniqResid( Atom *atom )
{
    memcpy(  atom->res_id, atom->chain, size_ID );
    sprintf( atom->res_id+size_ID, "%05d%c", 
                    atom->res_num_orig, atom->res_insert );
    memcpy(  atom->res_id+size_ID+6, atom->res_nam, size_ID );
    *(atom->res_id+2*size_ID+6) = '\0';
/*
    printf( "%s ", atom->res_id );
    PRINTID( atom->res_nam );
    putchar( '\n' );
*/
}

/****************  ReadUniqResid  ****************/
/*
**  
*/

void ReadUniqResid( Atom *atom )
{
    memcpy( atom->chain,   atom->res_id, size_ID );
    memcpy( atom->res_nam, atom->res_id+size_ID+6, size_ID );
    if ( sscanf(           atom->res_id+size_ID, "%5d%c", 
                 &atom->res_num_orig, &atom->res_insert ) < 2 )
    {
	printf( " * Cannot read from resid <%s>\n", atom->res_id );
	printf( " * Atom: " );
	WriteAtom( stdout, atom, use_both );
	putchar( '\n' );
	exit( 1 );
    }
}

/****************  WriteResidues  ****************/
/*
**  
*/

void WriteResidues( FILE *fp, Residue *residues, const int count )
{
    int    i = 0;

    while ( i++ < count )
    {
        WriteResidue( fp, residues );
        fprintf( fp, "# atoms: %2i", residues->atom_count );
        fputs( ". first at: ", fp );
        WriteAtom( fp, *(residues->atomptrs), use_ext );
        fputs( ". last at: ", fp );
        WriteAtom( fp, *((residues->atomptrs)+residues->atom_count-1), use_ext );
#ifdef DEBUG
	for ( j=0; j<residues->atom_count; j++ )
	    fprintf( fp, " %i", residues->atomidxs[j] );
#endif
        fputc( '\n', fp);
        residues++;
    }
}

/****************  CoordZero  ****************/
/*
**  
*/

void CoordZero( Coord *v )
{
    v->x = 0.F;
    v->y = 0.F;
    v->z = 0.F;
}

/****************  CoordAdd  ****************/
/*
**  
*/

void CoordAdd( Coord *v, const Coord c )
{
    v->x += c.x;
    v->y += c.y;
    v->z += c.z;
}

/****************  CoordDiv  ****************/
/*
**  
*/

void CoordDiv( Coord *v, const int n )
{
    v->x /= n;
    v->y /= n;
    v->z /= n;
}

/****************  FillAtom  ****************/
/*
**  Fill an atom with data. Only one of the two names is set.
**
**  NOTE: in principle the new, sequential residue number could be found
**  by going through the residue list, but for now comparisons will
**  be based on the residue ID (cf. FindAtom).
*/

Atom FillAtom( Resid res_id, ID res_nam, ID atom_nam, const int which )
{
/*    char   *fct_name = "FillAtom"; */
    Atom   at;

    InitAtom( &at );
    if ( which == use_int )
	memcpy( at.atom_int, atom_nam, size_ID );
    if ( which == use_ext )
	memcpy( at.atom_ext, atom_nam, size_ID );
    memcpy( at.res_nam, res_nam, size_ID );
    strcpy( at.res_id, res_id );
    ReadUniqResid( &at );
/*--at.res_num = at.res_num_orig;  */

/*
    STATUS("filling:");
    printf("%s ", res_id);
    PRINTIDv(res_nam);
    printf(" %i\n", at.res_num_orig);
*/
    return( at );
}

/****************  FindAtom1  ****************/
/*
**  Try to find an atom <at> in the external or internal (depending on <which>)
**  atom name list of a molecule <molptr>. 
**  Residues are matched on 'res_id'.
**  If <guess>, some name mappings are guessed (only for external names).
**  The list index is returned in <idx>, and the corresponding atom in
**  the function value.
*/

Atom *FindAtom1( int *idx, Atom *at, Mol *molptr, const int which, const Boolean guess )
{
/*    char     *fct_name = "FindAtom1"; */
    int      i = 0, count;
/*    int nmatch; */
    Atom     *atoms, *match, *try_variable;
    
    atoms = molptr->atoms;
    count = molptr->atom_count;

/* loop over all atoms and compare */

#ifdef DEBUG
    STATUS( "trying to find:" );
    printf( "   " );
    WriteAtom( stdout, at, use_both );
    putchar( '\n' );
#endif

    while ( i < count )
    {
#ifdef DEBUG
	printf( "looking at: " );
	WriteAtom( stdout, atoms, use_both );
#endif
/* match residue ID's */
        if ( EQ( atoms->res_id, at->res_id ) ) 
	{
/* match atom names */
            if ( which == use_ext )
	    {
		if ( EQID( atoms->atom_ext, at->atom_ext ) )
		    break;
	    }
            else if ( which == use_int )
	    {
		if ( EQID( atoms->atom_int, at->atom_int ) )
		    break;
	    }
	}
        atoms++;
        i++;
    }
    *idx = i;

/* found */

    if ( i < count )
        return( atoms );

/* not found */

    if ( which == use_ext && guess )
    {

/* derive a closely related name */

	if ( ( try_variable = GuessAtom( at ) ) == NULL )
	{
				    /* find wildcard matches */ 
#ifdef DEBUG
	    FindAtomsWC( &nmatch, at, molptr, use_ext );
#endif
	    return( NULL );
	}
	
/* try again */

#ifdef DEBUG
        STATUS( "trying again with:" );
	printf( "   " );
        WriteAtom( stdout, try_variable, use_both);
        putchar( '\n' );
#endif
	match = FindAtom1( idx, try_variable, molptr, use_ext, false );
	if ( match )
	{
	    printf( " External atom name remapped for: " );
	    WriteAtom( stdout, at, use_ext );
	    printf( " on: " );
	    PRINTIDv( try_variable->atom_ext );
	    putchar( '\n' );
	    return( match );
	}
	else
	{
				    /* find wildcard matches */ 
#ifdef DEBUG
	    FindAtomsWC( &nmatch, at, molptr, use_ext );
#endif
	    return( NULL );
	}
    }
    return( NULL );
}

/****************  FindAtom2  ****************/
/*
**  New version of FindAtom1.
**     - Faster: makes use of residue list.
**     - Correct residue names where needed.
**
**  Try to find an atom <at> in the external or internal (depending on <which>)
**  atom name list of a molecule <molptr>. 
**  Residues are matched on 'res_id'.
**  If <guess>, some name mappings are guessed (only for external names).
**  The list index is returned in <idx>, and the corresponding atom in
**  the function value.
**
**  If the residue list is scrambled  (i.e. atoms belonging to one residue are
**  scattered across more than one entry in molptr->residues) the routine 
**  fails (cf. comment in ExpandAtoms->MakeDefsMol).
*/

Atom *FindAtom2( int *idx, Atom *theAtom, Mol *molptr, const int which, 
                 const Boolean guess )
{
/*    char     *fct_name = "FindAtom2"; */
    Atom     *at = theAtom, *match, mod, *try_variable;
    Atomptr  *atomptrs;
    Residue  *res;
    char     *res_nam, *lastcharptr;
    int      *atomidxs, res_idx, count;
/*    int nmatch; */
    
#ifdef DEBUG
    STATUS( "trying to find:" );
    printf( "   " );
    WriteAtom( stdout, at, use_both );
    putchar( '\n' );
#endif

    find_atom_err = 0;

/* find the residue; map Arg+ onto Arg etc. */
/* 
   (The mapping is necessary because the expanded AtomLIB definitions will
   not contain Arg+ etc. (even though they exist in the unexpanded AtomLIB)
   if the set (i.e. the PDB file) does not contain Arg+ etc.; the mapping
   is only necessary to find the residue in the expanded AtomLIB.)
*/

    res_nam = at->res_nam;
    res_idx = FindResidueIndexFromId( at->res_id, molptr );
    if ( res_idx == -1 )
    {
	lastcharptr = LastNonBlank( at->res_nam, size_ID );
	if ( *lastcharptr == '+' || *lastcharptr == '-' )
	{
	    printf( " Residue name remapped for: " );
	    WriteAtom( stdout, at, which );
	    memcpy( &mod, at, sizeof( Atom ) );
	    at = &mod;
	    *LastNonBlank( at->res_nam, size_ID ) = ' ';
	    MakeUniqResid( at );
	    res_idx = FindResidueIndexFromId( at->res_id, molptr );
	    printf( " on: " );
	    PRINTIDv( at->res_nam );
	    putchar( '\n' );
	}
    }
    if ( res_idx == -1 )
    {
	find_atom_err = 1;
	return( NULL );
    }
#ifdef DEBUG
    sprintf( msg_string, "found residue index %i", res_idx );
    STATUSm;
#endif

/* find the atom */

    res = &molptr->residues[res_idx];
    count = res->atom_count;

    atomptrs = res->atomptrs;
    atomidxs = res->atomidxs;
    while ( count )
    {
#ifdef DEBUG
	printf( "looking at: " );
	WriteAtom( stdout, *atomptrs, use_both );
	printf( " index %i\n", *atomidxs );	    
#endif
	if ( which == use_ext )
	{
	    if ( EQID( (*atomptrs)->atom_ext, at->atom_ext ) )
	    {
		*idx = *atomidxs;
		break;
	    }
	}
	else if ( which == use_int )
	{
	    if ( EQID( (*atomptrs)->atom_int, at->atom_int ) )
	    {
		*idx = *atomidxs;
		break;
	    }
	}
        atomptrs++;
        atomidxs++;
	count--;
    }

/* found */

    if ( count )
        return( *atomptrs );

/* not found */

    if ( which == use_ext && guess )
    {

/* derive a closely related name */

	if ( ( try_variable = GuessAtom( at ) ) == NULL )
	{
				    /* find wildcard matches */ 
#ifdef DEBUG
	    FindAtomsWC( &nmatch, at, molptr, use_ext );
#endif
	    return( NULL );
	}
	
/* try again */

#ifdef DEBUG
        STATUS( "trying again with:" );
	printf( "   " );
        WriteAtom( stdout, try_variable, use_both);
        putchar( '\n' );
#endif
	match = FindAtom2( idx, try_variable, molptr, use_ext, false );
	if ( match )
	{
	    printf( " External atom name remapped for: " );
	    WriteAtom( stdout, at, use_ext );
	    printf( " on: " );
	    PRINTIDv( try_variable->atom_ext );
	    putchar( '\n' );
	    return( match );
	}
	else
	{
				    /* find wildcard matches */ 
#ifdef DEBUG
	    FindAtomsWC( &nmatch, at, molptr, use_ext );
#endif
	    return( NULL );
	}
    }
    return( NULL );
}

/****************  FindAtom3  ****************/
/*
**  New version of FindAtom1.
**  Handles chain and residue number wildcards such as set by AtomLIB.
**  This routine captures the matching of wildcards that was formerly
**  handled by ExpandAtoms, and is meant to be used by MapAtoms/MapPtrAtoms
**  only.
**  Residue name variants (such as Arg+) must be explicitly present in the
**  molecule stored at <molptr> (i.e. in the AtomLIB). 
**  In contrast to FindAtom2 the molecule does not have to contain a 
**  residue list. 
*/

Atom *FindAtom3( int *idx, Atom *at, Mol *molptr, const int which, 
                 const Boolean guess )
{
/*    char     *fct_name = "FindAtom3"; */
    int      i = 0, count;
/*    int nmatch; */
    Atom     *atoms, *match, *try_variable;
    
    atoms = molptr->atoms;
    count = molptr->atom_count;

/* loop over all atoms and compare */

#ifdef DEBUG
    STATUS( "trying to find:" );
    printf( "   " );
    WriteAtom( stdout, at, use_both );
    putchar( '\n' );
#endif

    while ( i < count )
    {
#ifdef DEBUG
	printf( "looking at: " );
	WriteAtom( stdout, atoms, use_both );
#endif
/* match residue and chain */
	if ( NOTEQID( atoms->res_nam, at->res_nam ) ) 
	{
	    atoms++;
	    i++;
	    continue;
	}
	if ( NOTEQID( atoms->chain, at->chain ) && atoms->chain[0] != '*' ) 
	{
	    atoms++;
	    i++;
	    continue;
	}
	if ( atoms->res_num_orig == at->res_num_orig ||
	     atoms->res_num_orig == wild_RESNUM )
	{
/* match atom names */
/* specific defs MUST come BEFORE wildcards in order to be found first */
            if ( which == use_ext )
	    {
		if ( EQID( atoms->atom_ext, at->atom_ext ) )
		    break;
	    }
            else if ( which == use_int )
	    {
		if ( EQID( atoms->atom_int, at->atom_int ) )
		    break;
	    }
	}
        atoms++;
        i++;
    }
    *idx = i;

/* found */

    if ( i < count )
        return( atoms );

/* not found */

    if ( which == use_ext && guess )
    {

/* derive a closely related name */

	if ( ( try_variable = GuessAtom( at ) ) == NULL )
	{
				    /* find wildcard matches */ 
#ifdef DEBUG
	    FindAtomsWC( &nmatch, at, molptr, use_ext );
#endif
	    return( NULL );
	}
	
/* try again */

#ifdef DEBUG
        STATUS( "trying again with:" );
	printf( "   " );
        WriteAtom( stdout, try_variable, use_both);
        putchar( '\n' );
#endif
	match = FindAtom3( idx, try_variable, molptr, use_ext, false );
	if ( match )
	{
	    printf( " External atom name remapped for: " );
	    WriteAtom( stdout, at, use_ext );
	    printf( " on: " );
	    PRINTIDv( try_variable->atom_ext );
	    putchar( '\n' );
	    return( match );
	}
	else
	{
				    /* find wildcard matches */ 
#ifdef DEBUG
	    FindAtomsWC( &nmatch, at, molptr, use_ext );
#endif
	    return( NULL );
	}
    }
    return( NULL );
}

/****************  TestFindAtom2  ****************/
/*
**  
*/

void TestFindAtom2( Mol *molptr )
{
    char     *fct_name = "TestFindAtom2";
    int      idx, count = molptr->atom_count;
    Atom     *atoms = molptr->atoms, *at;
    
    STATUS( "entering" );
    while ( count-- )
    {
	if ( ! ( at = FindAtom2( &idx, atoms, molptr, use_ext, true ) ) )
            WriteNotFoundAtom( stdout, atoms, use_ext, use_both );
	printf( "%i ", idx );
	if ( at != molptr->atoms + idx )
	{
	    puts( "error" );
	    WriteAtom( stdout, at, use_both ); putchar( '\n' );
	    WriteAtom( stdout, &atoms[idx], use_both ); putchar( '\n' );
	}
	else
	    putchar( '\n' );
	atoms++;
    }
    STATUS( "leaving" );
}

/****************  FindAtomsWC  ****************/
/*
**  Try to find an atom <at> in the external or internal (depending on <which>
**  atom name list of a molecule <molptr>, 
**  matching UP TO a wildcard character # or *. 
**  Residues are matched on 'res_id'.
**  The number of matches is returned in <nmatch>, and the address of an
**  Atomptr array in the function value.
*/

Atomptr *FindAtomsWC( int *nmatch, Atom *at, Mol *molptr, const int which )
{
/*    char     *fct_name = "FindAtomsWC"; */
    char     *wc;
    size_t   len = 0;
    char     test[size_ID+1];
    Atom     *atoms;
    Atomptr  *match = MatchList;
    int      i=0, n=0;
    int      count;
    
    atoms = molptr->atoms;
    count = molptr->atom_count;
    *nmatch = 0;

/* loop over all atoms and compare */

    if ( which == use_ext )
    {
	strncpy( test, at->atom_ext, size_ID );
	if ( ( wc = strpbrk( test, "#*" ) ) == NULL )
	    return( NULL );
	len = (size_t)( wc - test );
    }
    else if ( which == use_int )
    {
	strncpy( test, at->atom_int, size_ID );
	if ( ( wc = strpbrk( test, "#*" ) ) == NULL )
	    return( NULL );
	len = (size_t)( wc - test );
    }

    while ( i < count )
    {
/* match residue ID's -- assume unordered list */
        if ( EQ( atoms->res_id, at->res_id ) ) 
	{
/* match atom names */
            if ( which == use_ext )
	    {
		if ( EQn( atoms->atom_ext, at->atom_ext, len ) )
		{
		    if ( ++n > max_MATCH )
		    {
			puts( " * ERROR: max_MATCH too small" );
			exit( 1 );
		    }
		    *match++ = atoms;
		}
	    }
            else if ( which == use_int )
	    {
		if ( EQn( atoms->atom_int, at->atom_int, len ) )
		{
		    if ( ++n > max_MATCH )
		    {
			puts( " * ERROR: max_MATCH too small" );
			exit( 1 );
		    }
		    *match++ = atoms;
		}
	    }
	}
        atoms++;
        i++;
    }

    *nmatch = n;

#ifdef DEBUG
    if ( n ) {
	STATUS( "found wildcard representations:" );
	for ( i=0; i<n; i++ )
	{
	    printf( "   " );
	    WriteAtom( stdout, MatchList[i], use_both );
	    putchar('\n');
	}
    }
#endif

    if ( n  )
        return( MatchList );
    else
	return( NULL );
}

/****************  WriteNotFoundAtom  ****************/
/*
**  
*/

void WriteNotFoundAtom( FILE *fp, Atom *atom, const int find_which, 
                        const int print_which )
{
    if ( find_atom_err == 1 )
	fputs( " Residue name not found: ", fp );
    else
    {
	if ( find_which == use_ext )
	    fputs( " External atom name not found: ", fp );
	if ( find_which == use_int )
	    fputs( " Internal atom name not found: ", fp );
    }
    WriteAtom( fp, atom, print_which );
    fputc( '\n', fp );
}

/****************  FindResidueIndexFromNum  ****************/
/*
**  Returns the index in the residues-array of the residue 
**  for which the new residue number is <res_num>.
*/

int FindResidueIndexFromNum( const int res_num, Mol *molptr )
{
    int      i = 0, count;
    Residue  *residues;
    
    residues = molptr->residues;
    count = molptr->res_count;

    while ( i < count )
    {
	if ( res_num == residues->res_num_new )
	    return( i );
	residues++;
	i++;
    }
    return( -1 );
}

/****************  FindResidueIndexFromOrigNum  ****************/
/*
**  Returns the index in the residues-array of the residue 
**  for which the original residue number, chain id and insertion code are
**  <res_num>, <chain> and <insert>.
*/

int FindResidueIndexFromOrigNum( const int res_num, ID chain, char insert, 
                                 Mol *molptr )
{
    int      i = 0, count;
    Residue  *residues;
    
    residues = molptr->residues;
    count = molptr->res_count;

    while ( i < count )
    {
	if ( res_num == residues->res_num_orig &&
	     EQID( chain, residues->chain ) && insert == residues->res_insert )
	    return( i );
	residues++;
	i++;
    }
    return( -1 );
}

/****************  FindResidueIndexFromId  ****************/
/*
**  Returns the residue index of the residue 
**  for which the ID is <res_id>, by searching the residue list.
*/

int FindResidueIndexFromId( Resid res_id, Mol *molptr )
{
    int      i = 0, count;
    Residue  *residues;
    
    residues = molptr->residues;
    count = molptr->res_count;

    while ( i < count )
    {
	if ( EQ( res_id, residues->res_id ) )
	    return( i );
	residues++;
	i++;
    }
    return( -1 );
}

/****************  FindResidueNumFromId  ****************/
/*
**  Returns the new residue number of the residue 
**  for which the ID is <res_id>, by searching the residue list.
*/

int FindResidueNumFromId( Resid res_id, Mol *molptr )
{
    int      i = 0, count;
    Residue  *residues;
    
    residues = molptr->residues;
    count = molptr->res_count;

    while ( i < count )
    {
	if ( EQ( res_id, residues->res_id ) )
	    return( residues->res_num_new );
	residues++;
	i++;
    }
    return( -1 );
}

/****************  FindResidueNamFromId  ****************/
/*
**  Returns the residue name of the residue for which the ID
**  is <res_id> (not comparing the part which corresponds to the residue name 
**  of course...) by searching the residue list.
*/

char *FindResidueNamFromId( Resid res_id, Mol *molptr )
{
    int          i = 0, count;
    Residue      *residues;
    static char  unk[size_ID+1] = "UNK  ";
    
    residues = molptr->residues;
    count = molptr->res_count;

    while ( i < count )
    {
	if ( EQn( res_id, residues->res_id, size_ID+6 ) )
	    return( residues->res_nam );
	residues++;
	i++;
    }
    return( unk );
}

/****************  InitAtomS  ****************/
/*
**  
*/

void InitAtomS( Atom *atom, int count )
{
    while ( count-- )
        InitAtom( atom++ );
}

/****************  InitAtom  ****************/
/*
**  
*/

void InitAtom( Atom *atom )
{
    atom->res_num_new  = 0;
    atom->res_num_orig = 0;
    Blank( atom->chain,    size_ID );
    Blank( atom->atom_int, size_ID );
    Blank( atom->atom_ext, size_ID );
    Blank( atom->atom_element_symbol,  size_ID );
    Blank( atom->res_nam,  size_ID );
    atom->res_insert = ' ';
    *(atom->res_id) = '\0';
}

/****************  InitPsatomS  ****************/
/*
**  
*/

void InitPsatomS( Psatom *psatom, int count )
{
    while ( count-- )
        InitPsatom( psatom++ );
}

/****************  InitPsatom  ****************/
/*
**  
*/

void InitPsatom( Psatom *psatom )
{
    psatom->pseudo_typ   = 0;
    psatom->num_form     = 0;
    psatom->res_num_orig = 0;
    Blank( psatom->chain,      size_ID );
    Blank( psatom->res_nam,    size_ID );
    Blank( psatom->pseudo_nam, size_ID );
    Blank( (char *)psatom->form_nam,   size_ID*max_PSFORM );
    psatom->res_insert = ' ';
    *(psatom->res_id) = '\0';
}

/****************  GetLIBName  ****************/
/*
**  Determine the library name for the AtomLIB to be used.
**  If option <option> has been set (i.e. <option> is not item_NONE and 
**  the option value is not "none"), the routine uses the option value 
**  instead of presenting a menu.
*/

char *GetLIBName( const int option )
{
    Menu     *oldmenu;
    int      item;
    char     *name, *optval;

    static MenuItem menulist[] =
    {
        { item_aqua,   "aqua",      none, on }, 
        { item_disman, "disman",    none, on }, 
        { item_disgeo, "disgeo",    none, on }, 
        { item_diana,  "diana",     none, on }, 
        { item_xplor,  "xplor",     none, on }, 
        { item_xplori, "xplor-inv", none, on }, 
        { item_biosym, "biosym",    none, on }, 
        { item_whatif, "whatif",    none, on }, 
        { item_local,  "local",     none, on }, 
        { item_user,   "user",      none, on }, 
        { item_copy,   "copy",      none, on }, 
        { item_END,    "end",       none, on }  /* same as in AquaWhat, AquaPseudo and AquaHow */
    };
    static Menu menu =
    {
        "GetLIBName menu", "format", "Give atom name format:", 
        sizeof(menulist)/sizeof(MenuItem)
    };
    

    oldmenu = SetMenu( &menu, menulist );
    if ( option != item_NONE )
	optval = GetOptionCVal( option );       /* use value of option */
                                                /* OK if option menu has not been swapped */
    else
	optval = "none";

    while ( true )
    {
	if ( EQ( optval, "none" ) )
	    item = ReadMenuItem();
	else
	    item = FindMenuItem( optval );
        switch ( item )
        {
	    case item_aqua:
	    case item_disman:
	    case item_disgeo:
	    case item_diana:
	    case item_xplor:
	    case item_xplori:
	    case item_biosym:
	    case item_whatif:
	    case item_local:
	    case item_user:
	    case item_copy:
		name = GetMenuItemKey( item );
		ResetMenu( oldmenu );
		return( name );
            case item_END:
            case item_NULL:
                ResetMenu( oldmenu );
                return( NULL );	
            default:
                ProcessMenuError( menu_fatal );
                break;
	}
    }
}

/****************  ProcessAtomLIB  ****************/
/*
**  Read and process an AtomLIB i.e. a file with atom name mappings.
*/

Mol *ProcessAtomLIB( Mol *molptr, char *lib_nam, Filnam predef )
{
    char     *fct_name = "ProcessAtomLIB";
    char     *usrnam = NULL;
    Atom     *atmap_list;
    Mol      *atmap_mol;
    Filnam   lib_filnam;
    FILE     *lib_fp = NULL;
    Boolean  done = false;
    Menu     *oldmenu;
    int      atmap_count; 

    if ( EQ( lib_nam, "user" ) || EQ( lib_nam, "local" ) || EQ( lib_nam, "whatif" ) )

/* user defined library in Aqua format */
    {
	oldmenu = SetOneMenu( "Give name of atom-library file (or - to end):" );
	if ( EQ( lib_nam, "local" ) )
	    usrnam = "AtomLIB";
	else if ( EQ( lib_nam, "whatif" ) )
	    usrnam = predef;
	else
	    usrnam = GetOptionCVal( opt_atmlib );   /* use value of atomlib option */
                                                    /* OK if option menu has not been swapped */
	while ( !done )
	{
	    if ( !usrnam )
	    {
		MReadWord( lib_filnam );
		if ( *lib_filnam == '-' || strlen( lib_filnam ) == 0 )
		{
		    ResetMenu( oldmenu );
		    return( NULL );
		}
	    }
	    else
	    {
		strcpy( lib_filnam, usrnam );
		usrnam = NULL;                  /* to prevent infinite loop */
	    }
            if ( terseness < 6 ) {
	    printf( "Now opening: %s\n", lib_filnam );
            }
	    if ( ! ( done = OPENr( lib_fp, lib_filnam ) ) )
	    {
		printf( " * atom-library file <%s> could not be opened\n", lib_filnam );
		MError();
	    }
	}
	ResetMenu( oldmenu );
	if ( terseness < 6 ) {
	    STATUS( "reading user AtomLIB" );
        }
    }

    else

/* standard Aqua library from Aqua data directory */
    {
	strcpy( lib_filnam, "AtomLIB-" );
	strcat( lib_filnam, lib_nam );
	lib_fp = AquaDataFile( lib_filnam );
	if ( !lib_fp )
	    return( NULL );
        if ( terseness < 6 ) {
	    STATUS( "reading standard AtomLIB" );
        }
    }
    atmap_list = ReadAtomLIB( lib_fp, &atmap_count );
    if ( !atmap_list || atmap_count == 0 )
    {
	puts( " * Atom name definition list empty" );
	return( NULL );
    }

    sprintf( msg_string, "%i atom name definitions read", atmap_count );
    if ( terseness < 6 ) {
        STATUSm;
    }
    if ( OptionTrue( opt_prtlib ) )
	WriteFAtoms( stdout, atmap_list, atmap_count );


    if ( !OptionTrue( opt_find2 ) )
	return( MakeDefsMol( atmap_list, atmap_count, false ) );


/* expand the list from 'wildcard' residue numbers to individual residues */

    if ( terseness < 6 ) {
        STATUS( "generating definitions for present molecule" );
    }
    atmap_mol = ExpandAtoms( atmap_list, atmap_count, molptr );
    free( atmap_list );
    sprintf( msg_string, "%i atom name definitions active for present molecule", 
             atmap_mol->atom_count );
    if ( terseness < 6 ) {
        STATUSm;
    }

    if ( OptionTrue( opt_prtlib ) )
        WriteFAtoms( stdout, atmap_mol->atoms, atmap_mol->atom_count );

    return( atmap_mol );
}

/****************  ProcessStereoLIB  ****************/
/*
**  Read and process a StereoLIB i.e. a file describing the "equivalent
**  atoms" in prochiral groups.
**  This routine is an analogue of ProcessPseudo.
**  The descriptions and data types are the same as for pseudo atoms. 
*/

Psatom *ProcessStereoLIB( int *stereo_count, Mol *molptr, char *lib )
{
    char     *usrnam = NULL, *fct_name = "ProcessStereoLIB";
    char     *lib_nam;
    Psatom   *stereo_list;
    Filnam   lib_filnam;
    FILE     *lib_fp = NULL;
    Boolean  done = false;
    Menu     *oldmenu;

    lib_nam = GetLIBName( opt_rstinp );
    if ( !lib_nam )
	return( NULL );			         /* no processing done */

    if ( EQ( lib_nam, "user" ) )

/* type "user": user defined library in Aqua format */
    {
	oldmenu = SetOneMenu( "Give name of stereo-library file:" );
	usrnam = GetOptionCVal( opt_sterlib );   /* use value of stereolib option */
                                                 /* OK if option menu has not been swapped */
	while ( !done )
	{
	    if ( !usrnam )
		MReadWord( lib_filnam );
	    else
	    {
		strcpy( lib_filnam, usrnam );
		usrnam = NULL;                  /* to prevent infinite loop */
	    }

            if ( terseness < 6 ) {
	        printf( "Now opening: %s\n", lib_filnam );
            }
	    if ( ! ( done = OPENr( lib_fp, lib_filnam ) ) )
	    {
		printf( " * stereo-library file <%s> could not be opened\n", lib_filnam );
		MError();
	    }
	}
	ResetMenu( oldmenu );
        if ( terseness < 6 ) {
	    STATUS( "Reading user StereoLIB" );
        }
    }

    else

/* type "local", "diana", etc: standard Aqua library from Aqua data directory */
    {
	lib_fp = AquaDataFile( lib );
	if ( !lib_fp )
	    exit( 1 );
	    if ( terseness < 6 ) {
	        STATUS( "Reading standard StereoLIB" );
            }
    }
    stereo_list = ReadPseudoLIB( lib_fp, stereo_count );
    if ( …