distance constraint energy calculation in wham

[unres4.git] / source / wham / io_wham.F90

diff --git a/source/wham/io_wham.F90 b/source/wham/io_wham.F90
index 10364a4..187502e 100644 (file)
--- a/source/wham/io_wham.F90
+++ b/source/wham/io_wham.F90
@@ -461,11 +461,39 @@ allocate(ww(max_eneW))
       wscloc=ww(12)
       wtor=ww(13)
       wtor_d=ww(14)
+      wstrain=ww(15)
       wvdwpp=ww(16)
       wbond=ww(18)
       wsccor=ww(19)
       wcatcat=ww(42)
       wcatprot=ww(41)
+      wcorr3_nucl=ww(38)
+      wcorr_nucl=ww(37)
+      wtor_d_nucl=ww(36)
+      wtor_nucl=ww(35)
+      wsbloc=ww(34)
+      wang_nucl=ww(33)
+      wbond_nucl=ww(32)
+      welsb=ww(31)
+      wvdwsb=ww(30)
+      welpsb=ww(29)
+      wvdwpsb=ww(28)
+      welpp=ww(27)
+      wvdwpp_nucl=ww(26)
+      wscbase=ww(46)
+      wpepbase=ww(47)
+      wscpho=ww(48)
+      wpeppho=ww(49)
+!      print *,"KURWA",ww(48)
+!        "WSCBASE   ","WPEPBASE  ","WSCPHO    ","WPEPPHO   "
+!        "WVDWPP    ","WELPP     ","WVDWPSB   ","WELPSB    ","WVDWSB    ",&
+!        "WELSB     ","WBOND     ","WANG      ","WSBLOC    ","WTOR      ",&
+!        "WTORD     ","WCORR     ","WCORR3    ","WNULL     ","WNULL     ",&
+!        "WCATPROT  ","WCATCAT  
+!       +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+!       +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+!       +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+!       +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
 
       endif
 !
@@ -485,14 +513,32 @@ allocate(ww(max_eneW))
       weights(12)=wscloc
       weights(13)=wtor
       weights(14)=wtor_d
-      weights(15)=0 !wstrain !
-      weights(16)=0 !wvdwpp !
+      weights(15)=wstrain !0
+      weights(16)=wvdwpp !
       weights(17)=wbond
       weights(18)=0 !scal14 !
       weights(21)=wsccor
       weights(42)=wcatprot
       weights(41)=wcatcat
-
+      weights(26)=    wvdwpp_nucl 
+
+      weights(27) =welpp  
+      weights(28) =wvdwpsb
+      weights(29) =welpsb 
+      weights(30) =wvdwsb 
+      weights(31) =welsb  
+      weights(32) =wbond_nucl  
+      weights(33) =wang_nucl   
+      weights(34) =wsbloc 
+      weights(35) =wtor_nucl   
+      weights(36) =wtor_d_nucl 
+      weights(37) =wcorr_nucl  
+      weights(38) =wcorr3_nucl 
+      weights(41) =wcatcat
+      weights(42) =wcatprot
+      weights(46) =wscbase
+      weights(47) =wscpho
+      weights(48) =wpeppho
 ! el--------
       call card_concat(controlcard,.false.)
 
@@ -3682,8 +3728,9 @@ allocate(ww(max_eneW))
 !      include 'COMMON.HEADER'
 !      include 'COMMON.SBRIDGE'
       character(len=50) :: tytul
-      character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
-                      'D','E','F','G','H','I','J'/),shape(chainid))
+      character(len=1),dimension(24) :: chainid=reshape((/'A','B','C',&
+                      'D','E','F','G','H','I','J','K','L','M','N','O',&
+                      'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid))
       integer,dimension(nres) :: ica !(maxres)
       real(kind=8) :: temp,efree,etot,entropy,rmsdev
       integer :: ii,i,j,iti,ires,iatom,ichain,mnum
@@ -3700,9 +3747,11 @@ allocate(ww(max_eneW))
         mnum=molnum(i)
         iti=itype(i,mnum)
         if (iti.eq.ntyp1) then
+          if (itype(i-1,molnum(i-1)).eq.ntyp1) then
           ichain=ichain+1
           ires=0
           write (ipdb,'(a)') 'TER'
+          endif
         else
         ires=ires+1
         iatom=iatom+1