cmdmode.cc

c++编写的并行拉马克遗传算法的程序。实现分析对接程序

                        emap_total += emap[i];
                        elec_total += elec[i];
                        charge_total += charge[i];
                    } /*i*/
                    pr( logFile, "      __________  __________  __________  _______\n\n" );
                    pr( logFile, "Total %10.2f  %10.2f  %10.2f  %7.3f\n\n",
                        (emap_total + elec_total), emap_total, elec_total, charge_total );
                
                    pr( command_out_fp, "%.2f\n", etotal );
                    pr( logFile, "\n");
                    pr( logFile, "    E_estimated_free_energy_binding  = %.2f kcal/mol  [=(1)+(3)]\n", einter + torsFreeEnergy );
                    pr( logFile, "\n");
                    pr( logFile, "(1) E_intermolecular_non-bond        = %.2f kcal/mol\n", einter );
                    pr( logFile, "    E_intermolecular_atomic-affinity = %.2f kcal/mol\n", emap_total );
                    pr( logFile, "    E_intermolecular_electrostatic   = %.2f kcal/mol\n", elec_total );
                    pr( logFile, "(2) E_intramolecular_non-bond        = %.2f kcal/mol\n", eintra );
                    pr( logFile, "(3) E_intramolecular_torsional       = %.2f kcal/mol\n", torsFreeEnergy );
                    pr( logFile, "\n");
                    pr( logFile, "    E_docked                         = %.2f kcal/mol\n", einter+eintra );
                    pr( logFile, "\n");
                    fflush( logFile );
                    fflush( stderr );
                } /* endif */
                break;
/*
            ____________________________________________________________________
*/
            case COM_EVAL:
/*
                EVALUATE x y z  nx ny nz angle [in deg]
                tor1
                tor2
                ...
*/
                sscanf( command, "%*s %lf %lf %lf %lf %lf %lf %lf", 
                    &(S.T.x), &(S.T.y), &(S.T.z),  
                    &(S.Q.nx), &(S.Q.ny), &(S.Q.nz),  &(S.Q.ang) );
                pr( logFile, "COMMAND: eval %lf %lf %lf %lf %lf %lf %lf\n         ", 
                    S.T.x, S.T.y, S.T.z, S.Q.nx, S.Q.ny, S.Q.nz, S.Q.ang );
                /**/
                S.Q.ang = Rad( S.Q.ang );
                mkUnitQuat( &(S.Q) );
                /**/
                for (i=0; i<ntor; i++) {
                    fgets( command, LINE_LEN, command_in_fp );
                    sscanf( command, "%lf", &(S.tor[i]) ); /* S.tor in degrees */
                    pr( logFile, "%lf ", S.tor[i] ); /* S.tor in degrees */
                    S.tor[i] = Rad( S.tor[i] ); /* S.tor converted to radians */
                }
                cnv_state_to_coords( S,  vt, tlist, ntor,  crdpdb, crd, natom );
                if ( ntor > 0 ) {
                    eintra = eintcalPrint( nonbondlist, e_internal, crd, type, Nnb, B_calcIntElec, q1q2 );
                } else {
                    eintra = 0.0;
                }
                outside = FALSE;
                for (i = 0;  i < natom;  i++) {
                    outside = is_out_grid( crd[i][X], crd[i][Y], crd[i][Z] );
                    if ( outside ) {
                        break;
                    }
                } /*i*/
                if ( outside ) {
                    etotal = (einter = WallEnergy) + eintra;
                } else {
                    etotal = (einter = trilinterp( crd, charge, type, natom, map, inv_spacing, elec, emap, xlo, ylo, zlo )) + eintra;
                    /* lo ) ...using lo array is slower */
                }
                prStr( message, "%f\n", etotal );
                pr_2x( command_out_fp, logFile, message );
                fflush( logFile );
                fflush( stderr );
                break;
/*
            ____________________________________________________________________
*/
            case COM_OUTC:
                pr( logFile, "COMMAND: outc\n\n" );
                pr( logFile, "USER    Total Non-bonded  Energy of Complex = %+8.2f\n", etotal );
                for (i = 0;  i < natom;  i++) {
                    pr( command_out_fp, "%f %f %f\n",  crd[i][X], crd[i][Y], crd[i][Z] );
                    pr( logFile, "ATOM  %5d  %.8s%5d    %8.3f%8.3f%8.3f\n", i+1, rec8[i], nframes, crd[i][X], crd[i][Y], crd[i][Z] );
                } /*i*/
                strcpy( message, "TER\n\0" );
                pr_2x( command_out_fp, logFile, message );
                fflush( logFile );
                fflush( stderr );
                nframes++;
                break;
/*
            ____________________________________________________________________
*/
            case COM_OUTE:
                pr( logFile, "COMMAND: oute\n\n" );
                prStr( message, "USER    Total Internal Energy of Small Molecule = %.2f\n", eintra );
                pr_2x( command_out_fp, logFile, message );
                prStr( message, "USER    Total Docked Energy of Complex = %.2f\n", etotal );
                pr_2x( command_out_fp, logFile, message );
                prStr( message, "USER    Predicted Free Energy of Binding = %.2f\n", einter + torsFreeEnergy );
                pr_2x( command_out_fp, logFile, message );
                for (i = 0;  i < natom;  i++) {
                    pr( command_out_fp, "%.14s  %10.2f%10.2f\n", rec14[i], emap[i], elec[i] );
                    pr( logFile, "%.14s  %10.2f%10.2f\n", rec14[i], emap[i], elec[i] );
                } /*i*/
                fflush( logFile );
                fflush( stderr );
                break;
/*
            ____________________________________________________________________
*/
            case COM_TRJ:
/*
                Convert a trajectory, from .trj file into PDB...
*/
                sscanf( command, "%*s %s", trjFileName );
                pr( logFile, "COMMAND: traj %s\n\n", trjFileName );
                if ( !openFile( trjFileName, "r", &trjFile,jobStart,tms_jobStart,FALSE ) ) {
                    return -1;
                } else {
                    ntor_old = ntor;
                    T = 0.;

                    while( ( fgets( trjline, LINE_LEN, trjFile ) ) != NULL ) {
                        keyword_id = parse_trj_line( trjline );

                        if ( keyword_id == -1 ) {
                            pr( stderr, "%s: ERROR: Unrecognized keyword in Trajectory file: \"%s\"\n",
                                programname, trjline );
                            continue;
                        }

                        switch( keyword_id ) {
/*
                        ________________________________________________________
*/
                            case TRJ_NULL:
                                break;
/*
                        ________________________________________________________
*/
                            case TRJ_NTOR:
                                sscanf( trjline, "%*s %d", &ntor );
                                if ( ntor != ntor_old ) {
                                    prStr( message, "\n%s: mis-matched number of torsions in trajectory file %s - find correct docking-parameter file and substrate-PDBQ file.\n", programname, trjFileName );
                                    pr_2x( stderr, logFile, message );
                                }
                                break;
/*
                        ________________________________________________________
*/
                            case TRJ_RUN:
                                sscanf( trjline, "%*s %d", &irun );
                                break;
/*
                        ________________________________________________________
*/
                            case TRJ_CYCLE:
                                sscanf( trjline, "%*s %d", &icycle );
                                break;
/*
                        ________________________________________________________
*/
                            case TRJ_TEMP:
                                sscanf( trjline, "%*s %f", &T );
                                break;
/*
                        ________________________________________________________
*/
                            case TRJ_STATE:
                                if ( input_state( &S, trjFile, trjline, ntor, 
                                    &nstep, &E, &Eint, &lastmove ) != (int)0) {
                                    /*...input_state ensures tor is in radians*/
                                    cnv_state_to_coords( S, vt, tlist, ntor, 
                                        crdpdb, crd, natom);
                                    outside = FALSE;
                                    for (i = 0;  i < natom;  i++) {
                                        outside = is_out_grid( crd[i][X], 
                                            crd[i][Y], crd[i][Z] );
                                        if ( outside ) {
                                            break;
                                        }
                                    } /*i*/
                                    if ( outside ) {
                                        etotal = (einter = WallEnergy) + Eint;
                                    } else {
                                        etotal = (einter = trilinterp( crd, charge, type, natom, map, inv_spacing, elec, emap, xlo, ylo, zlo )) + Eint;
                                            /* lo ) + Eint; */
                                    } /* endif */
                                    pr( logFile, "USER   Run %d  Cycle %d  Step %d  Temp %.2f K  %c  Etot %.2f  Eint   %.2f\n", irun, icycle, nstep, T, lastmove, E, Eint   );
                                    switch ( lastmove ) {
                                        case 'A':
                                            movecode = 1;
                                            break;
                                        case 'a':
                                            movecode = 2;
                                            break;
                                        case 'e':
                                            movecode = 3;
                                            break;
                                        case 'R':
                                            movecode = 4;
                                            break;
                                        default:
                                            movecode = -1;
                                            break;
                                    } /* endswitch */
                                    for (i = 0;  i < natom;  i++) {
/*
Write coordinates and energy to
command output...

lastmove = AaeR,
movecode = 1234
*/
                                        pr( command_out_fp, "%.4f %.4f %.4f %.4f %.4f %.4f %.4f %.4f %.2f %d %d %d %d\n", crd[i][X], crd[i][Y], crd[i][Z],  emap[i], elec[i], emap[i]+elec[i],  Eint, etotal, T, movecode, nstep, icycle, irun );
                                        pr( logFile, "ATOM  %5d  %.8s%5d    %8.3f%8.3f%8.3f\n", i+1, rec8[i], nstep, crd[i][X], crd[i][Y], crd[i][Z] );
                                    } /*i*/
                                    pr( logFile, "TER\n" );
                                    fflush( logFile );
                                    fflush( stderr );
                                    nframes++;
                                } /* endif inputstate */
                                break;
/*
**                        ________________________________________________________
*/
                           default:
                                pr( stderr, "%s: ERROR: Unrecognized keyword: '%s'\n", programname, trjline );
                                break;
/*
**                        ________________________________________________________
*/
                        }

                    } /*  end of while -- EOF of trjFile  */
                    fclose( trjFile );
                }
                break;
/*
**            ____________________________________________________________________
*/
            default:
                pr( stderr, "%s: ERROR: Unrecognized command: '%s'\n", programname, line );
                break;
/*
**            ____________________________________________________________________
*/
        } /* End of switch. */
    } /* End of while. */
/*

Write out an AVS-readable field file, for input of trajectory file.

*/
    indpf = strindex( dock_param_fn, ".dpf" );
    strncpy( filename, dock_param_fn, ( size_t )indpf );
    filename[ indpf ] = '\0';
    strcat( filename, ".dlg.pdb\0" );

    for ( i=0, veclen = 13; i<veclen; i++ ) {
        offset[i] = i;
    }

    print_avsfld( logFile, veclen, natom, nframes, offset, 13, 
        "x y z E_atom E_elec E_atom_elec E_internal E_total Temp MoveCode Step Cycle Run", filename );

    jobEnd = times( &tms_jobEnd );
    timesys( jobEnd - jobStart, &tms_jobStart, &tms_jobEnd );
    pr_2x( logFile, stderr, UnderLine );

    return 0;
}
/* EOF */