|
#!/usr/bin/perl -w
#use strict;
#多基因的遗传算法
use AI::Genetic;
# benz_furan_MMTORSION
open FH,"benz_furan_mmenergyGA.tab" or die "can't find the file";
my @qmenergies=<FH>;
chomp(@qmenergies);
map($_=~s/^\d+\s+//,@qmenergies);
#
#map($_*=627.51,@qmenergies);
open FH,"benz_furan_mmenergyGA.tab";
my @mmenergies=<FH>;
chomp(@mmenergies);
map($_=~s/^\d+\s+//s,@mmenergies);
my @angles=(0 .. 72);
@angles=map($_*5,@angles);
my $ga = new AI::Genetic(
-fitness => \&fitnessFunc,
-type => 'rangevector',
-population => 1000,
-crossover => 0.85, #一个基因的话 没有交叉的必要0.9
-mutation => 0.1,
-terminate => \&terminateFunc,
);
$ga->init([
[1,300],
[1,8],
[0,36],
[0,300],
[1,8],
[0,36]
]);#基因的数目
$ga->evolve('rouletteSinglePoint', 1000); # 选择的方法 一共进化的代数
my $individual=$ga->getFittest(1);
my $resultx=$individual->genes();
my @resugenes=@$resultx;
print join "**",@resugenes;
print "\n";
my $resulty=$individual->score();
my $rmsdd=1/$resulty;
print " fittness: $resulty rmsd: $rmsdd\n";
sub fitnessFunc
{
my $genes = shift;
# print join " ", @$genes;
# die "111";
my @valuesmm=();#to save new mm energies
my $gene_pk1=@$genes[0]/10;
my $gene_pk2=@$genes[3]/10;
my $gene_pn1=@$genes[1];
my $gene_pn2=@$genes[4];
my $gene_phase1=@$genes[2]*10;
my $gene_phase2=@$genes[5]*10;
for(0..72)
{
# print "$mmenergies[$_]";
# print "$gene_pk1*(1+cos(($gene_pn1*$angles[$_]-$gene_phase1)/180*3.14159))\n";
# die "dd";
my $energydihe1=$gene_pk1*(1+cos(($gene_pn1*$angles[$_]-$gene_phase1)/180*3.14159));
my $energydihe2=$gene_pk2*(1+cos(($gene_pn2*$angles[$_]-$gene_phase2)/180*3.14159));
push @valuesmm,$energydihe2+$energydihe1+$mmenergies[$_];
}
my @valuesmm_sort=sort{$a <=> $b} @valuesmm;
@valuesmm=map($_-$valuesmm_sort[0],@valuesmm);
$rmsd=&calcrmsd(\@valuesmm,\@qmenergies);
#能用一个余弦函数表示就不用两个余弦函数表示;进行罚分处理
if($gene_pk1==$gene_pk2 and $gene_pn1 == $gene_pn1 and $gene_phase1 == $gene_phase2)
{
$rmsd*=2;
}
my $fitness;
$fitness=1/$rmsd;
print "$gene_pk1 $gene_phase1 $gene_pn1 $gene_pk2 $gene_pn1 $gene_phase2
fitness: $fitness $rmsd\n";
return $fitness;
}
sub calcrmsd
{
my @sander=@{$_[0]};
my @gauss=@{$_[1]};
my $rmsd=0;
if(@sander-0==@gauss-0)
{
for(0..@sander-1)
{
$rmsd+=($sander[$_]-$gauss[$_])**2;
}
return $rmsd/(@sander-0);
}
else
{
print "error ,the sander and qm have different energy\n";
exit 0;
}
}
sub terminateFunc
{
my $ga = shift;
# terminate if reached some threshold.
my $THRESHOLD=10000;
return 1 if $ga->getFittest->score > $THRESHOLD;
return 0;
}
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-11-1 07:14
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社