#!/usr/bin/perl
#
#-- Based on Backbone-enhanced ENM
#-- Generating DPA data for PDB protein structure
#     To predict Protein/Ligand interaction
#-- originally writen by DENGMING MING 
#     in Los Alamos National Lab, 2004
# 
#   ADD MULTIPLE-LEVEL SURFACE MSMS-POINTS ("-layer" option)
#   ADD REDUCED NMODES ("-nmodes" option)
#     by DENGMING MING
#     in Fudan University, 2013/11/29
#
#   ADD HIEARCHICAL CALC., I.E., GENERTING L_DPA-CLUSTERS FOR EACH (L)LEVEL,
#       AND FIND CONNECTION BETWEEN 1-, 2-,...,L-DPA-CLUSTERS
#       THE PREDICTION IS BASED ON THE CONNECTED-CLUSTERS 
#     in Nanjing Tech University, 2021/05/21
#
#   ADD MULTIPLE-CHAIN PROTEINS ("-chain" option)
#     in Nanjing Tech University, 2023/06/07
#                   
#
use strict;
use Sys::Hostname;
use lib "/xcommon/bin/_mdpa";
#use lib "/opt/xcommon/bin/_mdpa";

use _FILEprocc;
use _PDBprocc;
use _gDPA;
use _anaDPA;
use _gDCluster_DBindingsite;
use _cmpEXP;
#
#
#
my ($pdbid,$selechain,
    $cutcc,$delta_cutcc,$delta_cutlg,$cutlg,$cutll,$intcutcc,$wcc,$wctc,$wlg,$wll,$intwcc,
    $nmodes,$modekappa,$ev_threshold,$flag_enforceDPA,$flag_udpa,
    $msms_density,$msms_prob,$maxmsmslayer,$flag_layerconn,$connpoints,$conncutoff,$flag_recordlayer,
    $flag_msmsdense,$flag_msms_saveall,$flag_wmsms,$flag_wspdb,$flag_uMSMS,$input_msmsf,
    $epislon_dpa,$MinPts_dpa,$epislon_cluster_dpa,$cutdpaperct,$cutdpaperct2,$adjcutperct,
    $ndpavalue_section,$topdpa_num_cutoff,$dpabindingcutoff,$clu2ligcutoff,$flag_flexible_bcut,$flag_wcpdb,$flag_output_dpa,
    $flag_check,$flag_cmpexpri,$flag_wlbs,$flag_positive_pred,$flag_task,$flag_quite,
    $home,$scrh,$workdir,$exedir,$msmsdir,$structuredatadir,$outputdatadir,$dpadir,$totproccf
    )=&parse_arg();
#
print "\nMDPA_$pdbid> MULTIPLE-LAYERED DYNAMICS-PERTURBATION-ANALYSIS, NLAYER = $maxmsmslayer \n\n";
#
my $outputfile=$pdbid.'_processing.prd';
if($flag_recordlayer eq 'true'){
    open(MYOUTF,"> $outputdatadir/$outputfile") || die "CANNOT OPEN-WRITE DPADATA $outputdatadir/$outputfile\n";
}
my %mdpa_bindingsites=''; my %mdpa_centers=''; my %mdpa_nclu='';
undef %mdpa_bindingsites; undef %mdpa_centers; undef %mdpa_nclu;
my $nca=''; my $cacrd='';

my $flag_wspdb1=''; my $flag_wcpdb1='';
if($flag_recordlayer eq 'true'){
   $flag_wspdb1=$flag_wspdb; $flag_wcpdb1 = $flag_wcpdb;
}
my $msmslayer = 0;
while($msmslayer < $maxmsmslayer){
    $msmslayer++;
#    print "DPA> PROCESSING LAYER ----> $msmslayer :: $maxmsmslayer\n";
    my ($calc_info,$infocalc,$ndpaclu,$dpabindingsites,$dpacenters,
    $nca0,$cacrd0,$dpafile)=&_gDCluster_DBindingsite::protein_dpa($pdbid,$selechain,
    $cutcc,$delta_cutcc,$delta_cutlg,$cutlg,$cutll,$intcutcc,$wcc,$wctc,$wlg,$wll,$intwcc,
    $nmodes,$modekappa,$ev_threshold,$flag_enforceDPA,$flag_udpa,
    $msms_density,$msms_prob,$msmslayer,
    $flag_msmsdense,$flag_msms_saveall,$flag_wmsms,$flag_wspdb1,$flag_uMSMS,$input_msmsf,
    $epislon_dpa,$MinPts_dpa,$epislon_cluster_dpa,$cutdpaperct,$cutdpaperct2,$adjcutperct,
    $ndpavalue_section,$topdpa_num_cutoff,$dpabindingcutoff,$flag_flexible_bcut,$flag_wcpdb1,
    $flag_check,$flag_task,$flag_quite,
    $home,$scrh,$workdir,$exedir,$msmsdir,$structuredatadir,$outputdatadir,$dpadir,$totproccf
    );
    $mdpa_bindingsites{$msmslayer}=$dpabindingsites;
    $mdpa_centers{$msmslayer}=$dpacenters;
    $mdpa_nclu{$msmslayer}=$ndpaclu;
    $nca=$nca0; $cacrd=$cacrd0;
    if ($flag_output_dpa eq 'false'){ system("rm  $outputdatadir"."/"."$dpafile")}
    else {print "DPA_CALC_$pdbid> DPA recorded: $outputdatadir/$dpafile\n" if $flag_check eq 'true';}
    if($flag_recordlayer eq 'true'){
        foreach my $rkid (sort {$a<=>$b} keys %$dpabindingsites){
#            print "000SITE> $msmslayer -- $rkid --> @{$$dpabindingsites{$rkid}}\n";
            foreach my $crdid (sort {$a<=>$b} @{$$dpabindingsites{$rkid}}){
                my $resName=@{$$cacrd{$crdid}}[3]; my $chainID=@{$$cacrd{$crdid}}[4];
                my $resSeq =@{$$cacrd{$crdid}}[5]; my $iCode  =@{$$cacrd{$crdid}}[6];
    #           printf  ">>>> %3d %3d %3s %1s %5d %1s\n",$msmslayer,$rkid,$resName,$chainID,$resSeq,$iCode;
                printf MYOUTF "%3d --> %3d :: %3s %1s %5d %1s\n",$msmslayer,$rkid,$resName,$chainID,$resSeq,$iCode;
            }
        }
    }
}
close(MYOUTF) if($flag_recordlayer eq 'true');
 
#
#-- CONECTIVITY-ANALYSIS OF CLUSTERS-AT-DIFFERENT-LAYERS


my ($nclu,$predictions,$clusters,$conn,$noise_centers)=&_anaDPA::get_layercluster_connection($pdbid,$maxmsmslayer,
          $connpoints,$conncutoff,$epislon_dpa,$MinPts_dpa,$epislon_cluster_dpa,
          \%mdpa_nclu,\%mdpa_bindingsites,\%mdpa_centers,$nca,$cacrd,$flag_check);

#   foreach my $rkid (sort {$a<=>$b} keys %{$predictions}){
#        print "DDD_CONN> $pdbid --> $rkid ---> $$conn{$rkid}\n";
#        print "DDD_SITE> $pdbid --> $rkid ---> @{$$predictions{$rkid}}\n";
#   }
#   print "CENTERS>>>>>>>>>>>>>>> \n";
#   foreach my $rkid (sort {$a<=>$b} keys %{$clusters}){
#        foreach my $crdid (sort {$a<=>$b} keys %{$$clusters{$rkid}}){
#            print "DDD_CTRS> $pdbid --> $rkid -- $crdid: @{$$clusters{$rkid}{$crdid}}\n"; 
#        }
#   }

my $outputfile=$pdbid.'.prd';
open(MYOUTF1,"> $outputdatadir/$outputfile") || 
	die "CANNOT OPEN-WRITE PREDICTIONS $outputdatadir/$outputfile\n";
    foreach my $rkid (sort {$a<=>$b} keys %{$conn}){
        printf MYOUTF1 "#%3d :: %s\n",$rkid,$$conn{$rkid};}
    foreach my $rkid (sort {$a<=>$b} keys %$predictions){
        foreach my $crdid (sort {$a<=>$b} @{$$predictions{$rkid}}){
            my $resName=@{$$cacrd{$crdid}}[3]; my $chainID=@{$$cacrd{$crdid}}[4];
            my $resSeq =@{$$cacrd{$crdid}}[5]; my $iCode  =@{$$cacrd{$crdid}}[6];
    #        printf  ">>>> %3d %3d %3s %1s %5d %1s\n",$maxmsmslayer,$rkid,$resName,$chainID,$resSeq,$iCode;
            printf MYOUTF1 " %3d :: %3s %1s %5d %1s\n",$rkid,$resName,$chainID,$resSeq,$iCode;
    }}
close(MYOUTF1);
if (-e "$outputdatadir/$outputfile"){
    print "\nMDPA_$pdbid> DPA-BINDINGSITE RECORDED: $outputdatadir/$outputfile\n" if $flag_check eq 'true';
}

if($flag_wcpdb eq 'true'){
    my $outputfile=$pdbid.'_dpaclusters.pdb';
    open(MYOUTF2,"> $outputdatadir/$outputfile") || 
        die "CANNOT OPEN-WRITE DPACLUSTERS $outputdatadir/$outputfile\n";
    foreach my $rkid (sort {$a<=>$b} keys %{$conn}){
        printf MYOUTF2 "REMARK  %3d :: %s\n",$rkid,$$conn{$rkid};}
    my $iatom=0; #my @chainid_list=[O,P,Q,R,S,T,U,V,W,X,Y,Z]; 
    my $chainid='';my @chainid_list=("O".."Z"); my $resname='';
    foreach my $cid (sort {$a<=>$b} keys %$clusters) {
        $chainid=@chainid_list[($cid-1)%scalar(@chainid_list)] if($cid >=1);
        $resname='TOP';$resname='NOS' if($cid eq 0);
        $chainid = "x" if($cid eq 0);
        foreach my $i (sort {$a<=>$b} keys %{$$clusters{$cid}}){
            $iatom++; my  $resseq=$i;
            my $x = @{$$clusters{$cid}{$i}}[0];
            my $y = @{$$clusters{$cid}{$i}}[1];
            my $z = @{$$clusters{$cid}{$i}}[2];
            my $pid =@{$$clusters{$cid}{$i}}[-1]; $pid =9999 if $pid >= 9999;

           #printf  "%6s%5d %4s%1s%3s %1s%4d%1s   %8.3f%8.3f%8.3f%6.2f%6.2f%6s%4s%4d\n", 'HETATM',
           #$iatom,' CA ',' ',$resname,$chainid,$resseq,' ',$x,$y,$z,1.0,10.0,'','MING',$pid;
           printf MYOUTF2 "%6s%5d %4s%1s%3s %1s%4d%1s   %8.3f%8.3f%8.3f%6.2f%6.2f%6s%4s%4d\n", 'HETATM',
           $iatom,' CA ',' ',$resname,$chainid,$resseq,' ',$x,$y,$z,1.0,10.0,'','MING',$pid;
            }
        }
    close(MYOUTF2);
    if (-e "$outputdatadir/$outputfile"){
        print "MDPA_$pdbid> DPA-CLUSTER-PDB RECORDED: $outputdatadir/$outputfile\n";
    }
}

exit() if $flag_cmpexpri ne 'true';
&_cmpEXP::cmp_expri($pdbid,$selechain,$structuredatadir,$workdir,$clu2ligcutoff,
			$nca,$conn,$predictions,$clusters,$cacrd,$maxmsmslayer,
			$flag_wlbs,$flag_positive_pred,$flag_check);

exit();

sub parse_arg{
    my ($pdbid,$pdbf);
    my ($cutcc,$delta_cutcc,$delta_cutlg,$cutlg,$cutll,$intcutcc,$wcc,$wctc,$wlg,$wll,$intwcc);
    my ($nmodes,$modekappa,$ev_threshold,$flag_enforceDPA,$flag_udpa,$flag_quickdpa);
    my ($maxmsmslayer,$flag_layerconn,$connpoints,$conncutoff,$flag_recordlayer);
    my ($msms_density,$msms_prob,$flag_msmsdense,$flag_msms_saveall,$flag_wmsms,$flag_wspdb,$flag_uMSMS,$input_msmsf);
    my ($flag_check,$flag_cmpexpri,$flag_wlbs,$flag_positive_pred);
    my ($epislon_dpa,$MinPts_dpa,$epislon_cluster_dpa);
    my ($cutdpaperct,$cutdpaperct2,$adjcutperct,$ndpavalue_section,$topdpa_num_cutoff,$dpabindingcutoff,$clu2ligcutoff,$flag_flexible_bcut,$flag_wcpdb,$flag_output_dpa);
    my ($home,$scrh,$workdir,$exedir,$msmsdir,$structuredatadir,$outputdatadir,$dpadir,$totproccf);

    my $totproccf='_DPAproccessing.log';

    my $numArgs = $#ARGV + 1; #&show_usage if $numArgs <= 0;
    my $index=-1;
    while (++$index<$numArgs){
	my $para=$ARGV[$index];
	if($para eq "-p") {
	    $pdbid=$ARGV[++$index];} 
	elsif($para eq "-f") {
	    $pdbf=$ARGV[++$index];} 

	elsif($para eq "-chain"){
	    do {
		$selechain=$selechain.$ARGV[$index]." " if( $ARGV[++$index] !~ /^-/);
	    } until($index>=$numArgs || $ARGV[$index] =~ /^-/);
	    $index-=1 if( $ARGV[$index] =~ /^-/);
	}

#	elsif($para eq "-dpa"){           #PARAMETERS FOR GENERATING DPA VALUES
#	    $dpastart=$ARGV[++$index]; 
#	    $dpastep=$ARGV[++$index]; }

#
#-- dynamics
	elsif($para eq "-wcc"){
	    $wcc=$ARGV[++$index];}
	elsif($para eq "-wctc"){
	    $wctc=$ARGV[++$index];}
	elsif($para eq "-wlg"){
	    $wlg=$ARGV[++$index];}
	elsif($para eq "-wll"){
	    $wll=$ARGV[++$index];}
	elsif($para eq "-intwcc"){
	    $intwcc=$ARGV[++$index];}
	elsif($para eq "-cutcc"){
	    $cutcc=$ARGV[++$index];}
	elsif($para eq "-dcutcc"){
	    $delta_cutcc=$ARGV[++$index];}
	elsif($para eq "-cutlg"){
	    $cutlg=$ARGV[++$index];}
	elsif($para eq "-cutll"){
	    $cutll=$ARGV[++$index];}
	elsif($para eq "-dcutlg"){
	    $delta_cutlg=$ARGV[++$index];}
	elsif($para eq "-intcutcc"){
	    $intcutcc=$ARGV[++$index];}
	elsif($para eq "-threshold" or $para eq "-thres"){
	    $ev_threshold=$ARGV[++$index];}
	elsif($para eq "-nmodes"){
	    $nmodes=$ARGV[++$index];}	
#	elsif($para eq "-modekappa"){
#	    $modekappa=$ARGV[++$index];}
	elsif($para eq "-quickdpa"){
	    $flag_quickdpa='true';}	

	elsif($para eq "-enforcedpa" or $para eq "-edpa"){
	    $flag_enforceDPA='true';}      #also write out CLUSTER info in PDB
	elsif($para eq "-rundpa" or $para eq "-rdpa"){
	    $flag_udpa='false';}
#
#--MSMS surface poitns
	elsif($para eq "-msms"){             #parameter to run MSMS
	    $msms_density=$ARGV[++$index];
	    $msms_prob=$ARGV[++$index];}
	elsif($para eq "-msmslayer" or $para eq "-layer"){ #define maximum layers of MSMS-points on protein surface
	    $maxmsmslayer=$ARGV[++$index];}                   #layer 2 is built on layer 1, and so on
	elsif($para eq "-layerconn"){      #analysis the connection of DPA-clusters in different layers 
            $flag_layerconn='true';}	   #new feature added by Dengming Ming, 2021.5. NJTECH.
        elsif($para eq "-connpoints" or $para eq "-npconn"){
            $connpoints=$ARGV[++$index];}  
        elsif($para eq "-cutoffconn"){
            $conncutoff=$ARGV[++$index];}  

	elsif($para eq "-msmsdense" or $para eq "-dense"){
	    $flag_msmsdense='true';}

	elsif($para eq "-msmsdense2" or $para eq "-dense2"){
	    $flag_msmsdense=$ARGV[++$index];}

	elsif($para eq "-msmssaveall"){              #parameter to run MSMS
	    $flag_msms_saveall='true';}

	elsif($para eq "-umsms"){             #reading MSMS points from external file
	    $flag_uMSMS='true';
	    $input_msmsf=$ARGV[++$index]; }

	elsif($para eq "-showlayer"){
	    $flag_recordlayer='true';}

	elsif($para eq "-ws" or $para eq "-wmsms"){
	    $flag_wmsms='true';}

	elsif($para eq "-wspdb"){
	    $flag_wspdb='true';}
#
#--DPA clusters & predictions
	elsif($para eq "-clup"){             #OPTICS CLUSTERING PARAMETERS
	    $epislon_dpa=$ARGV[++$index];
	    $MinPts_dpa=$ARGV[++$index]; 
	    $epislon_cluster_dpa=$ARGV[++$index];} 

	elsif($para eq "-topp"){
	    $cutdpaperct=$ARGV[++$index];}  #TOP percent threshold in EVD fitting	
	elsif($para eq "-topp2"){
	    $cutdpaperct2=$ARGV[++$index];}  #direct TOP percent threshold
	elsif($para eq "-adjcutperct" and $para eq "-perct"){
	    $adjcutperct='true';}             #Make sure at least  one NON-noise cluster is generated by changing $cutdpaperct
	elsif($para eq "-topnum" or $para eq "-topn"){
	    $topdpa_num_cutoff=$ARGV[++$index];}
	elsif($para eq "-binding_cutoff" or $para eq "-bcut"){
	    $dpabindingcutoff=$ARGV[++$index];}
	elsif($para eq "-flexible-bcut" or $para eq "-fbcut"){
            $flag_flexible_bcut='true'}
	elsif($para eq "-clear-dpaf"){
	    $flag_output_dpa='false';}
	elsif($para eq "-wcpdb"){
	    $flag_wcpdb='true';}          #also write out CLUSTER info in PDB
	elsif($para eq "-clu2ligcutoff" or $para eq "-ccut"){
	    $clu2ligcutoff=$ARGV[++$index];}

	elsif($para eq "-clearp"){
	    $flag_positive_pred='true';}          #also write out CLUSTER info in PDB

#
#-- controls
	elsif($para eq "-check"){
	    $flag_check='true';}

	elsif($para eq "-cmpexpri"){
	    $flag_cmpexpri='true';}

	elsif($para eq "-wlbs"){
	    $flag_wlbs='true';}

	elsif($para eq "-task"){
	    $flag_task=$ARGV[++$index];
	    $flag_udpa='false' if ($flag_task eq 'genedpa' or $flag_task eq 'gdpa');}	

	elsif($para eq "-quite" or $para eq "-q" ){
	    $flag_quite='true';}
#
##--others
	elsif($para eq "-h") {
	    &show_usage;}

	else{
	    print "\nUnrecognizable parameter: \"$para\", \"-h\" for help \n\n";
	    &processing_record($outputdatadir,'error_inp',$totproccf,-1,"\nbad parameter for \"$para\", \"-h\" for help ");
	    exit();}
    }

#
#
#-- input files/data
    $selechain='all' if(! $selechain);   #in default DPA use the whole protein
#
##-- dynamics parameters
##--
    $cutcc=10 if ! $cutcc;
    $delta_cutcc=1. if ! $delta_cutcc;  #ignore this parameter ?
    $delta_cutlg=5. if ! $delta_cutlg;
    $cutlg=$cutcc+$delta_cutlg if ! $cutlg or  $cutlg < $cutcc+$delta_cutlg;  
    $cutll=13. if ! $cutll;

    $wcc=1. if ! $wcc;
    $wctc=40. if ! $wctc;
    $wlg=12.  if ! $wlg ;
    $wll=10.  if ! $wll ;

    $intwcc=$wcc if (! $intwcc);
    $intcutcc=$cutcc if(! $intcutcc);

    $nmodes=-1 if ! $nmodes;   #number of modes for DPA calc., "-1" for 3N-modes
    $nmodes=-9 if $flag_quickdpa eq 'true';
    
    $modekappa=100 if ! $modekappa;         #egeinvector kappa value in selecting modes
    $ev_threshold=0.001 if ! $ev_threshold; #maximum positive value for zero-modes eigenvalues
    $flag_enforceDPA='false' if ! $flag_enforceDPA;  #enforce DPA calc. on weird structures
    $flag_udpa='true' if ! $flag_udpa;    #(true/false) Use/NOT_use pre-calculated DPA data

#
##-- MSMS parameters
#
    $msms_density=1. if ! $msms_density;  #NO INPUT FOR DENSITY IN DEFAULT
    $msms_prob=1.4   if ! $msms_prob;     #NO INPUT FOR PROB IN DEFAULT
    $maxmsmslayer=3     if ! $maxmsmslayer; #generate first layer MSMS points
    if ($maxmsmslayer <= 0){
        print "DPA_WARNING> msmslayer was set as 1, ignore zero/negative input: $maxmsmslayer\n";
        $maxmsmslayer = 1}
    $flag_msms_saveall='false' if ! $flag_msms_saveall;       #do not save MSMS points of intermediate layers 
    $flag_wmsms='false' if ! $flag_wmsms;  #do not write MSMS points in MSMS format
    $flag_wspdb='false' if ! $flag_wspdb;  #do not write MSMS points in PDB format
    $input_msmsf='' if ! $input_msmsf;     #EXTERNAL MSMS FILE AS INPUT
    if (! $flag_uMSMS) {              #DO NOT USE EXTERNAL MSMS FILE
	$flag_uMSMS='false';
	$input_msmsf='';
    }
    $flag_msmsdense='false' if ! $flag_msmsdense;
#
#--DPA CLUSTERS

#-  Clusterization
    $epislon_dpa=5 if ! $epislon_dpa;  #Clusterization: OPTICS parameters
    $MinPts_dpa=3 if ! $MinPts_dpa;
    $epislon_cluster_dpa=5 if ! $epislon_cluster_dpa;    
    $cutdpaperct=0.96 if ! $cutdpaperct;      #threshold to selet TOP DPA points in EVD fitting
    $cutdpaperct2=0.99 if ! $cutdpaperct2;     #threshold used in selet TOP DPA points
    $adjcutperct='false' if ! $adjcutperct;    #Make sure at least  one NON-noise cluster is generated by decreasing $cutdpaperct
    $ndpavalue_section=40 if ! $ndpavalue_section; #section number used in dividing DPA values to make profile density 
    $topdpa_num_cutoff=5 if ! $topdpa_num_cutoff;  #the threshold in selecting TOP_DPA clusters
    $dpabindingcutoff=6. if ! $dpabindingcutoff;  #cutoff in searching residues from DPA clusters
    $clu2ligcutoff=3.5  if ! $clu2ligcutoff;    #cutoff to search ligand in a pocket
    $flag_flexible_bcut='false' if ! $flag_flexible_bcut;  #DO NOT change bcut for different LAYER-clusters
    $flag_wcpdb='false' if ! $flag_wcpdb;        #DO NOT WRITE DPA CLUSTERS IN PDB FORMAT
    $flag_check='false' if ! $flag_check;
    $flag_cmpexpri='false' if ! $flag_cmpexpri;

#
#-- layer connectivity
    $flag_layerconn='true' if ! $flag_layerconn;  #do layer connectivity analysis
    $connpoints=-1 if ! $connpoints;              # connecitvy between clu1/2 determination:
                                                  # -1 using OPTICS to determin 
                                                  # n (>=1) check if exists n-pairwise contacts.
    if ($connpoints ne "-1" && $connpoints !~ /^(\d+)$/){
        print "\nPlease enter \"-1\" or positive integer for -npconn ($connpoints), -h for help\n\n"; exit;}
    $conncutoff=$epislon_cluster_dpa if ! $conncutoff;
    
#
#-- Directories
    my $host = hostname; $host= uc($host); 
    $home='.';     $scrh='.';     $workdir=$scrh;

    $exedir='/xcommon/bin/_mdpa/exe'; 
    $msmsdir='/xcommon/lib/msms/2.6.1';
#    $exedir='/opt/xcommon/bin/_mdpa/exe'; 
#    $msmsdir='/opt/xcommon/bin/_mdpa/exe/msms/2.6.1';

    $structuredatadir='.';
    $outputdatadir='.';
    $dpadir='.';

    if($pdbid =~ /(.*)\.pdb/){ $pdbid = $1}
    if($pdbf =~ /(.*)\.pdb/){ $pdbid = $1}
    if (! $pdbid and ! $pdbf ){
        print "\nPlease input protein-file using option \"-p\" or \"-f\", -h for help\n\n";exit;}

    return $pdbid,$selechain,
    $cutcc,$delta_cutcc,$delta_cutlg,$cutlg,$cutll,$intcutcc,$wcc,$wctc,$wlg,$wll,$intwcc,
    $nmodes,$modekappa,$ev_threshold,$flag_enforceDPA,$flag_udpa,
    $msms_density,$msms_prob,$maxmsmslayer,$flag_layerconn,$connpoints,$conncutoff,$flag_recordlayer,
    $flag_msmsdense,$flag_msms_saveall,$flag_wmsms,$flag_wspdb,$flag_uMSMS,$input_msmsf,
    $epislon_dpa,$MinPts_dpa,$epislon_cluster_dpa,$cutdpaperct,$cutdpaperct2,$adjcutperct,
    $ndpavalue_section,$topdpa_num_cutoff,$dpabindingcutoff,$clu2ligcutoff,$flag_flexible_bcut,$flag_wcpdb,$flag_output_dpa,
    $flag_check,$flag_cmpexpri,$flag_wlbs,$flag_positive_pred,$flag_task,$flag_quite,
    $home,$scrh,$workdir,$exedir,$msmsdir,$structuredatadir,$outputdatadir,$dpadir,$totproccf
}

sub show_usage{

    print "\nDynamic Perturbation Method for Predition Functional Sites On Protein Surface\n\n";
    print "USAGE: dpa -protein-option  [-dynamics-option] [-msms-option] [-cluster-option] [-calc_dpa-option]\n";
    print "                            [-directory-option]  [-task-option] [-control-option]\n\n";
    print "  -protein-option:      -p  pdbfilename xxxx.pdb (also for -f) or its headname\n";
    print "                    -chain  selected chainid (default all-chains)\n";
    print "\n";
    print " -dynamics-option:    -wcc  interaction strength (gamma=1.0) between C-alpha's NOT immediately jointed\n";  
    print "                     -wctc  interaction strength (gamma=40.0) between two C-alpha atoms immediately jointed\n"; 
    print "                      -wlg  interaction strength (gamma=12.0) between ligand and C-alpha atoms\n"; 
    print "                      -wll  interaction strength (gamma=10.0) among ligand atoms\n";  
    print "                   -intwcc  interface (chain-chain) gamma  (=wcc)\n";
    print "                    -cutcc  interaction cutoff (10.0) among C-alpha atoms\n"; 
    print "                    -cutlg  interaction cutoff (15.0) between ligand and C-alpha atoms\n"; 
    print "                    -cutll  interaction cutoff (13.0) among ligand atoms\n"; 
    print "                 -intcutcc  interface cutoff (=cutcc)\n";
    print "                    -thres  maximum positive eignvalue (0.001) allowed for zero-modes\n"; 
    print "                   -nmodes  # of modes in DPA calc. (-1 for 3N-modes),(-5: N/2; -9: N, -99: 2N)\n";
#    print "                -modekappa  negeinvector kappa value in selecting mode (100)\n";
    print "                 -quickdpa  flag for quick DPA calc. by setting nmodes: -9 for N-modes\n";
    print "\n";
    print "     -msms-option:   -msms  density,prob-radius for running MSMS program\n";
    print "                    -umsms  msms file in MSMS format, using msms file\n";
    print "                -msmslayer  (=1) maximum MSMS-point layer, layer 2 built on 1, and so on (-layer)\n";
    print "                -layerconn  analyze connections between DPA-clusters in different layers\n";
    print "               -connpoints  (=-1, also -npconn) mini # of pairwised point-contacts at \"connected\" clusters\n";
    print "                            -1 :: connectivity determined by OPTICS alg. using \"-clup parameters \"\n";
    print "               -cutoffconn  (10) maximum distance defining two contacted points\n";
    print "                -showlayer  record DPA-Site predictions for each layer clusters\n";
    print "                    -dense  generate most dense-distribution surface points\n";
    print "                   -dense2  genearte denser point sets (double/triple/...)\n";
    print "              -msmssaveall  recording MSMS points of all layers: xxxx_i.msms xxxx_msms_i.pdb\n";
    print "                     -surf  output file for MSMS points format: x,y,z,dx,dy,dz,id\n";
    print "                    -wspdb  output file for MSMS points in PDB-format\n";
    print "\n";  
    print "  -cluster-option:  -clup   OPTICS CLUSTERING parameters: epislon(5), minimum_points(3), episolon_cluster(5)\n";
    print "\n";
    print "-directory-option:\n";
    print "\n";  
    print " -calc_dpa-option:   -topp  percentage threshold for TOP MSMS points (0.96 =96%), used in EVD fitting\n";
    print "                    -topp2  lowest percentage allowed in searching TOP MSMS points(99%)\n";
    print "                    -perct  or -adjcutperct:  threshold of TOP percent (reduced) in EVD fitting\n";
    print "                     -edpa  or -enforcedpa:  enforce DPA running even EVD fitting failed with -topp value\n";
    print "                     -rdpa  or -rundpa:  enforce running DPA and ignored existed dpa-file\n";
    print "                    -wcpdb  write out DPA cluster points in PDB format\n";
    print "               -clear-dpaf  keep DPA data files (false)\n";
    print "                     -topn  minimum number of TOP POINTS set for DPA prediction ANALYSIS\n";
    print "                     -bcut  cutoff (6.0) to search function sites from TOP MSMS points (or -binding_cutoff)\n";
    print "                     -ccut  cutoff (3.5) to search ligand from TOP MSMS points (or -clu2ligcutoff)\n";
    print "                    -fbcut  increase bcut :  bcut =bcut0 + msmsprob*(layer -1 )\n";
    print "\n";
    print "     -task-option:  -task   \"genedpa/gdpa\" for generate DPA data ONLY, else for normal DPA prediction calc.\n";
    print "\n";
    print "  -control-option:  -check  reserve intemediate files for check\n";
    print "                 -cmpexpri  perform comparison with experiment LBS (PDB-SITE data)\n";
    print "                     -wlbs  write PDB SITE info (like \"gsite\"\n";
    print "                   -clearp  ignore clusters with negative MCC prediction\n";
    print "                    -quite  suppress screen print (or -q)\n";
    print "\n";
    exit;
}