remarks file: nmr_relaxation/refine_new.inp
remarks Relaxation matrix refinement using the method of
remarks  Yip and Case (1989) 

set seed 584930 end                  {*Seed for random velocity assignment.*}

structure @c2.psf end                                {*Read structure file.*}
coor @c2_01_sa.pdb                                      {*Read coordinates.*}

parameter                                          {*Set energy parameters.*}
    @TOPPAR:parallhdg.pro
        nbonds                 
        atom  nbxmod -3  tolerance 0.5  
        repel 0.75  rcon 4.0  rexp 2  cutnb 4.5   
    end
end
flags include relaxation end          {*Include the relaxation energy term.*}

relaxation
    nrestraints 200                        {*Expected number of restraints.*}

{*========================== Specify the physical parameters of the system.*}

    group      DMSO                 {*A group of spectra, DMSO, is defined.*}
    unresolved DMSO                      
        methyl                             {*Methyl protons are unresolved.*}
    end
    average 3                               {*R-3 average to methyl groups.*}
    omega      DMSO 500.0e6                {*Spectrometer frequency, in Hz.*}
    taucorrel  DMSO                              {*Correlation times, in s.*}
        model lipari                                    {*Order parameters.*}
        vector (name hn or name ha) 
               (name hn or name ha)                     0.75e-9 0.85 
        vector (name hn or name ha) 
               (name h* and not (name hn or name ha))   0.75e-9 0.80 
        vector (name h* and not (name hn or name ha)) 
               (name h* and not (name hn or name ha))   0.75e-9 0.65
    end 

    classification D300                    {*Define a classification, D300,*}
    @c2_300.int 

    clgroup     D300 DMSO                  {*which belongs to group DMSO.  *}
    taumix      D300 0.30                              {*Mixing time, in s.*}
    minint      D300 0.05                   {*Minimum measurable intensity.*}

    cutoff                {*Cutoff used for gradient and relaxation matrix.*}
        mode all  value 4.5
    end

    potential parabola                      {*Error estimates are not used.*}
    iexp 0.1666666667                     {*Dev = I_calc^iexp - I_obs^iexp.*}
    eexp 2                                                  {*E = Dev^eexp.*}

    calibrate         {*Calculate spectrum and turn automatic rescaling on.*}
        quality 0.33  automatic on  reference all  group class  
    end

    print  threshold 9999  end        {*Print R value of initial structure.*} 
end

{*============================================== Cool from 1000 K to 100 K.*}

vector do (mass = 10) (all)                               {*Uniform masses.*}
vector do (fbeta= 50) (all)   {*Friction coefficient for heatbath coupling.*}

evaluate ($init_t = 1500)            {*Initial annealing temperature, in K.*}
evaluate ($final_t =  75)              {*Final annealing temperature, in K.*}
evaluate ($tempstep = 25)               {*Temperature step per cycle, in K.*}

vector do (vx = maxwell($init_t)) (all)        {*Assign initial velocities.*}
vector do (vy = maxwell($init_t)) (all)
vector do (vz = maxwell($init_t)) (all)

evaluate ($totalstep = 2000)              {*Total annealing time, in steps.*}
evaluate ($nstep = int($totalstep*$tempstep/($init_t-$final_t)))

relaxation 
    weight    D300  1000                       {*Set the relaxation weight.*}
    tolerance 0.0              {*Calculate relaxation energy at every step.*} 
end

evaluate ($bath = $init_t)
while ($bath > $final_t) loop temp
    dynamics  verlet
        nstep=$nstep timest=0.001 {ps}
        iasvel=current  
        tcoupling=true  tbath=$bath
        nprint=$nstep iprfrq=$nstep
    end
    evaluate ($bath = $bath - $tempstep)     {*Reduce heatbath temperature.*}
end loop temp

{*======================================== Conjugate gradient minimization.*}

evaluate ($count = 0)
while ($count < 3) loop mini
    evaluate ($count = $count + 1)
    relaxation
        calibrate                    {*Automatic calibration has to be off.*}
           quality 0.33  automatic off  reference all  group class  
        end
    end   
    minimize powell nstep 30 end 
end loop mini

write coor output= refine.pdb end 

stop