zmiana do DiL
[unres.git] / source / unres / src_CSA_DiL / parmread.F
1       subroutine parmread
2 C
3 C Read the parameters of the probability distributions of the virtual-bond
4 C valence angles and the side chains and energy parameters.
5 C
6 C Important! Energy-term weights ARE NOT read here; they are read from the
7 C main input file instead, because NO defaults have yet been set for these
8 C parameters.
9 C
10       implicit real*8 (a-h,o-z)
11       include 'DIMENSIONS'
12 #ifdef MPI
13       include "mpif.h"
14       integer IERROR
15 #endif
16       include 'COMMON.IOUNITS'
17       include 'COMMON.CHAIN'
18       include 'COMMON.INTERACT'
19       include 'COMMON.GEO'
20       include 'COMMON.LOCAL'
21       include 'COMMON.TORSION'
22       include 'COMMON.SCCOR'
23       include 'COMMON.SCROT'
24       include 'COMMON.FFIELD'
25       include 'COMMON.NAMES'
26       include 'COMMON.SBRIDGE'
27       include 'COMMON.MD_'
28       include 'COMMON.SETUP'
29       character*1 t1,t2,t3
30       character*1 onelett(4) /"G","A","P","D"/
31       character*1 toronelet(-2:2) /"p","a","G","A","P"/
32       logical lprint,LaTeX
33       dimension blower(3,3,maxlob)
34       dimension b(13)
35       character*3 lancuch,ucase
36 C
37 C For printing parameters after they are read set the following in the UNRES
38 C C-shell script:
39 C
40 C setenv PRINT_PARM YES
41 C
42 C To print parameters in LaTeX format rather than as ASCII tables:
43 C
44 C setenv LATEX YES
45 C
46       call getenv_loc("PRINT_PARM",lancuch)
47       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
48       call getenv_loc("LATEX",lancuch)
49       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
50 C
51       dwa16=2.0d0**(1.0d0/6.0d0)
52       itypro=20
53 C Assign virtual-bond length
54       vbl=3.8D0
55       vblinv=1.0D0/vbl
56       vblinv2=vblinv*vblinv
57 c
58 c Read the virtual-bond parameters, masses, and moments of inertia
59 c and Stokes' radii of the peptide group and side chains
60 c
61 #ifdef CRYST_BOND
62       read (ibond,*) vbldp0,akp,mp,ip,pstok
63       do i=1,ntyp
64         nbondterm(i)=1
65         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
66         dsc(i) = vbldsc0(1,i)
67         if (i.eq.10) then
68           dsc_inv(i)=0.0D0
69         else
70           dsc_inv(i)=1.0D0/dsc(i)
71         endif
72       enddo
73 #else
74       read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
75       do i=1,ntyp
76        read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
77      &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
78         dsc(i) = vbldsc0(1,i)
79         if (i.eq.10) then
80           dsc_inv(i)=0.0D0
81         else
82           dsc_inv(i)=1.0D0/dsc(i)
83         endif
84       enddo
85 #endif
86       if (lprint) then
87         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
88         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
89      &   'inertia','Pstok'
90         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
91         do i=1,ntyp
92           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
93      &      vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
94           do j=2,nbondterm(i)
95             write (iout,'(13x,3f10.5)')
96      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
97           enddo
98         enddo
99       endif
100 #ifdef CRYST_THETA
101 C
102 C Read the parameters of the probability distribution/energy expression 
103 C of the virtual-bond valence angles theta
104 C
105       do i=1,ntyp
106         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
107      &    (bthet(j,i),j=1,2)
108         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
109         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
110         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
111         sigc0(i)=sigc0(i)**2
112       enddo
113       close (ithep)
114       if (lprint) then
115       if (.not.LaTeX) then
116         write (iout,'(a)') 
117      &    'Parameters of the virtual-bond valence angles:'
118         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
119      & '    ATHETA0   ','         A1   ','        A2    ',
120      & '        B1    ','         B2   '        
121         do i=1,ntyp
122           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
123      &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
124         enddo
125         write (iout,'(/a/9x,5a/79(1h-))') 
126      & 'Parameters of the expression for sigma(theta_c):',
127      & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
128      & '     ALPH3    ','    SIGMA0C   '        
129         do i=1,ntyp
130           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
131      &      (polthet(j,i),j=0,3),sigc0(i) 
132         enddo
133         write (iout,'(/a/9x,5a/79(1h-))') 
134      & 'Parameters of the second gaussian:',
135      & '    THETA0    ','     SIGMA0   ','        G1    ',
136      & '        G2    ','         G3   '        
137         do i=1,ntyp
138           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
139      &       sig0(i),(gthet(j,i),j=1,3)
140         enddo
141        else
142         write (iout,'(a)') 
143      &    'Parameters of the virtual-bond valence angles:'
144         write (iout,'(/a/9x,5a/79(1h-))') 
145      & 'Coefficients of expansion',
146      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
147      & '   b1*10^1    ','    b2*10^1   '        
148         do i=1,ntyp
149           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
150      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
151         enddo
152         write (iout,'(/a/9x,5a/79(1h-))') 
153      & 'Parameters of the expression for sigma(theta_c):',
154      & ' alpha0       ','  alph1       ',' alph2        ',
155      & ' alhp3        ','   sigma0c    '        
156         do i=1,ntyp
157           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
158      &      (polthet(j,i),j=0,3),sigc0(i) 
159         enddo
160         write (iout,'(/a/9x,5a/79(1h-))') 
161      & 'Parameters of the second gaussian:',
162      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
163      & '        G2    ','   G3*10^1    '        
164         do i=1,ntyp
165           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
166      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
167         enddo
168       endif
169       endif
170 #else 
171
172 C Read the parameters of Utheta determined from ab initio surfaces
173 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
174 C
175       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
176      &  ntheterm3,nsingle,ndouble
177       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
178       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
179       do i=1,maxthetyp
180         do j=1,maxthetyp
181           do k=1,maxthetyp
182             aa0thet(i,j,k)=0.0d0
183             do l=1,ntheterm
184               aathet(l,i,j,k)=0.0d0
185             enddo
186             do l=1,ntheterm2
187               do m=1,nsingle
188                 bbthet(m,l,i,j,k)=0.0d0
189                 ccthet(m,l,i,j,k)=0.0d0
190                 ddthet(m,l,i,j,k)=0.0d0
191                 eethet(m,l,i,j,k)=0.0d0
192               enddo
193             enddo
194             do l=1,ntheterm3
195               do m=1,ndouble
196                 do mm=1,ndouble
197                  ffthet(mm,m,l,i,j,k)=0.0d0
198                  ggthet(mm,m,l,i,j,k)=0.0d0
199                 enddo
200               enddo
201             enddo
202           enddo
203         enddo
204       enddo 
205       do i=1,nthetyp
206         do j=1,nthetyp
207           do k=1,nthetyp
208             read (ithep,'(3a)',end=111,err=111) res1,res2,res3
209             read (ithep,*,end=111,err=111) aa0thet(i,j,k)
210             read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
211             read (ithep,*,end=111,err=111)
212      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
213      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
214      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
215      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
216             read (ithep,*,end=111,err=111)
217      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
218      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
219      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
220           enddo
221         enddo
222       enddo
223 C
224 C For dummy ends assign glycine-type coefficients of theta-only terms; the
225 C coefficients of theta-and-gamma-dependent terms are zero.
226 C
227       do i=1,nthetyp
228         do j=1,nthetyp
229           do l=1,ntheterm
230             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
231             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
232           enddo
233           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
234           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
235         enddo
236         do l=1,ntheterm
237           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
238         enddo
239         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
240       enddo
241 C
242 C Control printout of the coefficients of virtual-bond-angle potentials
243 C
244       if (lprint) then
245         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
246         do i=1,nthetyp+1
247           do j=1,nthetyp+1
248             do k=1,nthetyp+1
249               write (iout,'(//4a)') 
250      &         'Type ',onelett(i),onelett(j),onelett(k) 
251               write (iout,'(//a,10x,a)') " l","a[l]"
252               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
253               write (iout,'(i2,1pe15.5)')
254      &           (l,aathet(l,i,j,k),l=1,ntheterm)
255             do l=1,ntheterm2
256               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') 
257      &          "b",l,"c",l,"d",l,"e",l
258               do m=1,nsingle
259                 write (iout,'(i2,4(1pe15.5))') m,
260      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
261      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
262               enddo
263             enddo
264             do l=1,ntheterm3
265               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
266      &          "f+",l,"f-",l,"g+",l,"g-",l
267               do m=2,ndouble
268                 do n=1,m-1
269                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
270      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
271      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
272                 enddo
273               enddo
274             enddo
275           enddo
276         enddo
277       enddo
278       call flush(iout)
279       endif
280       write (2,*) "Start reading THETA_PDB"
281       do i=1,ntyp
282         read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
283      &    (bthet(j,i),j=1,2)
284         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
285         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
286         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
287         sigc0(i)=sigc0(i)**2
288       write (2,*) "End reading THETA_PDB"
289       enddo
290       close (ithep_pdb)
291 #endif
292       close(ithep)
293 #ifdef CRYST_SC
294 C
295 C Read the parameters of the probability distribution/energy expression
296 C of the side chains.
297 C
298       do i=1,ntyp
299         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
300         if (i.eq.10) then
301           dsc_inv(i)=0.0D0
302         else
303           dsc_inv(i)=1.0D0/dsc(i)
304         endif
305         if (i.ne.10) then
306         do j=1,nlob(i)
307           do k=1,3
308             do l=1,3
309               blower(l,k,j)=0.0D0
310             enddo
311           enddo
312         enddo  
313         bsc(1,i)=0.0D0
314         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
315      &    ((blower(k,l,1),l=1,k),k=1,3)
316         do j=2,nlob(i)
317           read (irotam,*,end=112,err=112) bsc(j,i)
318           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
319      &                                 ((blower(k,l,j),l=1,k),k=1,3)
320         enddo
321         do j=1,nlob(i)
322           do k=1,3
323             do l=1,k
324               akl=0.0D0
325               do m=1,3
326                 akl=akl+blower(k,m,j)*blower(l,m,j)
327               enddo
328               gaussc(k,l,j,i)=akl
329               gaussc(l,k,j,i)=akl
330             enddo
331           enddo 
332         enddo
333         endif
334       enddo
335       close (irotam)
336       if (lprint) then
337         write (iout,'(/a)') 'Parameters of side-chain local geometry'
338         do i=1,ntyp
339           nlobi=nlob(i)
340           if (nlobi.gt.0) then
341             if (LaTeX) then
342               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
343      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
344                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
345      &                             'log h',(bsc(j,i),j=1,nlobi)
346                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
347      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
348               do k=1,3
349                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
350      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
351               enddo
352             else
353               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
354               write (iout,'(a,f10.4,4(16x,f10.4))')
355      &                             'Center  ',(bsc(j,i),j=1,nlobi)
356               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
357      &           j=1,nlobi)
358               write (iout,'(a)')
359             endif
360           endif
361         enddo
362       endif
363 #else
364
365 C Read scrot parameters for potentials determined from all-atom AM1 calculations
366 C added by Urszula Kozlowska 07/11/2007
367 C
368       do i=1,ntyp
369         read (irotam,*,end=112,err=112) 
370        if (i.eq.10) then 
371          read (irotam,*,end=112,err=112) 
372        else
373          do j=1,65
374            read(irotam,*,end=112,err=112) sc_parmin(j,i)
375          enddo  
376        endif
377       enddo
378 C
379 C Read the parameters of the probability distribution/energy expression
380 C of the side chains.
381 C
382       do i=1,ntyp
383         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
384         if (i.eq.10) then
385           dsc_inv(i)=0.0D0
386         else
387           dsc_inv(i)=1.0D0/dsc(i)
388         endif
389         if (i.ne.10) then
390         do j=1,nlob(i)
391           do k=1,3
392             do l=1,3
393               blower(l,k,j)=0.0D0
394             enddo
395           enddo
396         enddo  
397         bsc(1,i)=0.0D0
398         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
399      &    ((blower(k,l,1),l=1,k),k=1,3)
400         do j=2,nlob(i)
401           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
402           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
403      &                                 ((blower(k,l,j),l=1,k),k=1,3)
404         enddo
405         do j=1,nlob(i)
406           do k=1,3
407             do l=1,k
408               akl=0.0D0
409               do m=1,3
410                 akl=akl+blower(k,m,j)*blower(l,m,j)
411               enddo
412               gaussc(k,l,j,i)=akl
413               gaussc(l,k,j,i)=akl
414             enddo
415           enddo 
416         enddo
417         endif
418       enddo
419       close (irotam_pdb)
420 #endif
421       close(irotam)
422
423 #ifdef CRYST_TOR
424 C
425 C Read torsional parameters in old format
426 C
427       read (itorp,*,end=113,err=113) ntortyp,nterm_old
428       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
429       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
430       do i=1,ntortyp
431         do j=1,ntortyp
432           read (itorp,'(a)')
433           do k=1,nterm_old
434             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i) 
435           enddo
436         enddo
437       enddo
438       close (itorp)
439       if (lprint) then
440         write (iout,'(/a/)') 'Torsional constants:'
441         do i=1,ntortyp
442           do j=1,ntortyp
443             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
444             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
445           enddo
446         enddo
447       endif
448 #else
449 C
450 C Read torsional parameters
451 C
452       read (itorp,*,end=113,err=113) ntortyp
453       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
454       do iblock=1,2
455       do i=-ntyp,-1
456        itortyp(i)=-itortyp(-1)
457       enddo
458 c      write (iout,*) 'ntortyp',ntortyp
459       do i=0,ntortyp,ntortyp-1
460         do j=-ntortyp,ntortyp
461           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
462      &          nlor(i,j,iblock)
463            nterm(-i,-j,iblock)=nterm(i,j,iblock)
464            nlor(-i,-j,iblock)=nlor(i,j,iblock)
465           v0ij=0.0d0
466           si=-1.0d0
467           do k=1,nterm(i,j,iblock)
468             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
469      &      v2(k,i,j,iblock)
470             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
471             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock) 
472             v0ij=v0ij+si*v1(k,i,j,iblock)
473             si=-si
474           enddo
475           do k=1,nlor(i,j,iblock)
476             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
477      &        vlor2(k,i,j),vlor3(k,i,j) 
478             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
479           enddo
480           v0(i,j)=v0ij
481         enddo
482       enddo
483       close (itorp)
484       if (lprint) then
485         write (iout,'(/a/)') 'Torsional constants:'
486         do i=1,ntortyp
487           do j=1,ntortyp
488             write (iout,*) 'ityp',i,' jtyp',j
489             write (iout,*) 'Fourier constants'
490             do k=1,nterm(i,j,iblock)
491               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
492      &         v2(k,i,j,iblock)
493             enddo
494             write (iout,*) 'Lorenz constants'
495             do k=1,nlor(i,j,iblock)
496               write (iout,'(3(1pe15.5))') 
497      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
498             enddo
499           enddo
500         enddo
501       endif
502 C
503 C 6/23/01 Read parameters for double torsionals
504 C
505       do iblock=1,2
506       do i=0,ntortyp-1
507         do j=-ntortyp+1,ntortyp-1
508           do k=-ntortyp+1,ntortyp-1
509             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
510               write (iout,*) "OK onelett",
511      &         i,j,k,t1,t2,t3
512
513             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) 
514      &        .or. t3.ne.toronelet(k)) then
515               write (iout,*) "Error in double torsional parameter file",
516      &         i,j,k,t1,t2,t3
517 #ifdef MPI
518               call MPI_Finalize(Ierror)
519 #endif
520                stop "Error in double torsional parameter file"
521             endif
522             read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
523      &         ntermd_2(i,j,k,iblock)
524             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
525      &         ntermd_1(i,j,k,iblock))
526             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
527      &         ntermd_1(i,j,k,iblock))
528             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
529      &         ntermd_1(i,j,k,iblock))
530             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
531      &         ntermd_1(i,j,k,iblock))
532 C Matrix of D parameters for one dimesional foureir series
533             do l=1,  ntermd_1(i,j,k,iblock)
534              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
535              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
536              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
537              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
538             enddo
539             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
540      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
541      &         v2s(m,l,i,j,k,iblock),
542      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
543 C Matrix of D parameters for two dimesional fourier series
544             do l=1,ntermd_2(i,j,k,iblock)
545              do m=1,l-1
546              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
547              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
548              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
549              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
550              enddo!m
551             enddo!l
552           enddo!k
553         enddo!j
554       enddo!i
555       enddo!iblock
556 cc      if (lprint) then
557       write (iout,*) 
558       write (iout,*) 'Constants for double torsionals'
559       do i=1,ntortyp
560         do j=-ntortyp,ntortyp
561           do k=-ntortyp,ntortyp
562             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
563      &        ' nsingle',ntermd_1(i,j,k,iblock),
564      &        ' ndouble',ntermd_2(i,j,k,iblock)
565             write (iout,*)
566             write (iout,*) 'Single angles:'
567             do l=1,ntermd_1(i,j,k,iblock)
568               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
569      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
570      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
571             enddo
572             write (iout,*)
573             write (iout,*) 'Pairs of angles:'
574             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
575             do l=1,ntermd_2(i,j,k,iblock)
576               write (iout,'(i5,20f10.5)') 
577      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
578             enddo
579             write (iout,*)
580             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
581             do l=1,ntermd_2(i,j,k,iblock)
582               write (iout,'(i5,20f10.5)') 
583      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
584      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
585             enddo
586             write (iout,*)
587           enddo
588         enddo
589       enddo
590 cc      endif
591 #endif
592 C
593 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
594 C         correlation energies.
595 C
596       read (isccor,*,end=119,err=119) nterm_sccor
597       do i=1,20
598         do j=1,20
599           read (isccor,'(a)')
600           do k=1,nterm_sccor
601             read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
602      &        v2sccor(k,i,j) 
603           enddo
604         enddo
605       enddo
606       close (isccor)
607       if (lprint) then
608         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
609         do i=1,20
610           do j=1,20
611             write (iout,*) 'ityp',i,' jtyp',j
612             do k=1,nterm_sccor
613               write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
614             enddo
615           enddo
616         enddo
617       endif
618 C
619 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
620 C         interaction energy of the Gly, Ala, and Pro prototypes.
621 C
622       if (lprint) then
623         write (iout,*)
624         write (iout,*) "Coefficients of the cumulants"
625       endif
626       read (ifourier,*) nloctyp
627       do i=1,nloctyp
628       do i=0,nloctyp-1
629         read (ifourier,*,end=115,err=115)
630         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
631         if (lprint) then
632         write (iout,*) 'Type',i
633         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
634         endif
635         B1(1,i)  = b(3)
636         B1(2,i)  = b(5)
637 c        b1(1,i)=0.0d0
638 c        b1(2,i)=0.0d0
639         B1tilde(1,i) = b(3)
640         B1tilde(2,i) =-b(5) 
641 c        b1tilde(1,i)=0.0d0
642 c        b1tilde(2,i)=0.0d0
643         B2(1,i)  = b(2)
644         B2(2,i)  = b(4)
645 c        b2(1,i)=0.0d0
646 c        b2(2,i)=0.0d0
647         CC(1,1,i)= b(7)
648         CC(2,2,i)=-b(7)
649         CC(2,1,i)= b(9)
650         CC(1,2,i)= b(9)
651 c        CC(1,1,i)=0.0d0
652 c        CC(2,2,i)=0.0d0
653 c        CC(2,1,i)=0.0d0
654 c        CC(1,2,i)=0.0d0
655         Ctilde(1,1,i)=b(7)
656         Ctilde(1,2,i)=b(9)
657         Ctilde(2,1,i)=-b(9)
658         Ctilde(2,2,i)=b(7)
659 c        Ctilde(1,1,i)=0.0d0
660 c        Ctilde(1,2,i)=0.0d0
661 c        Ctilde(2,1,i)=0.0d0
662 c        Ctilde(2,2,i)=0.0d0
663         DD(1,1,i)= b(6)
664         DD(2,2,i)=-b(6)
665         DD(2,1,i)= b(8)
666         DD(1,2,i)= b(8)
667 c        DD(1,1,i)=0.0d0
668 c        DD(2,2,i)=0.0d0
669 c        DD(2,1,i)=0.0d0
670 c        DD(1,2,i)=0.0d0
671         Dtilde(1,1,i)=b(6)
672         Dtilde(1,2,i)=b(8)
673         Dtilde(2,1,i)=-b(8)
674         Dtilde(2,2,i)=b(6)
675 c        Dtilde(1,1,i)=0.0d0
676 c        Dtilde(1,2,i)=0.0d0
677 c        Dtilde(2,1,i)=0.0d0
678 c        Dtilde(2,2,i)=0.0d0
679         EE(1,1,i)= b(10)+b(11)
680         EE(2,2,i)=-b(10)+b(11)
681         EE(2,1,i)= b(12)-b(13)
682         EE(1,2,i)= b(12)+b(13)
683 c        ee(1,1,i)=1.0d0
684 c        ee(2,2,i)=1.0d0
685 c        ee(2,1,i)=0.0d0
686 c        ee(1,2,i)=0.0d0
687 c        ee(2,1,i)=ee(1,2,i)
688       enddo
689       if (lprint) then
690       do i=1,nloctyp
691         write (iout,*) 'Type',i
692         write (iout,*) 'B1'
693         write(iout,*) B1(1,i),B1(2,i)
694         write (iout,*) 'B2'
695         write(iout,*) B2(1,i),B2(2,i)
696         write (iout,*) 'CC'
697         do j=1,2
698           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
699         enddo
700         write(iout,*) 'DD'
701         do j=1,2
702           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
703         enddo
704         write(iout,*) 'EE'
705         do j=1,2
706           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
707         enddo
708       enddo
709       endif
710
711 C Read electrostatic-interaction parameters
712 C
713       if (lprint) then
714         write (iout,*)
715         write (iout,'(/a)') 'Electrostatic interaction constants:'
716         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
717      &            'IT','JT','APP','BPP','AEL6','AEL3'
718       endif
719       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
720       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
721       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
722       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
723       close (ielep)
724       do i=1,2
725         do j=1,2
726         rri=rpp(i,j)**6
727         app (i,j)=epp(i,j)*rri*rri 
728         bpp (i,j)=-2.0D0*epp(i,j)*rri
729         ael6(i,j)=elpp6(i,j)*4.2D0**6
730         ael3(i,j)=elpp3(i,j)*4.2D0**3
731         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
732      &                    ael6(i,j),ael3(i,j)
733         enddo
734       enddo
735 C
736 C Read side-chain interaction parameters.
737 C
738       read (isidep,*,end=117,err=117) ipot,expon
739       if (ipot.lt.1 .or. ipot.gt.5) then
740         write (iout,'(2a)') 'Error while reading SC interaction',
741      &               'potential file - unknown potential type.'
742 #ifdef MPI
743         call MPI_Finalize(Ierror)
744 #endif
745         stop
746       endif
747       expon2=expon/2
748       if(me.eq.king)
749      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
750      & ', exponents are ',expon,2*expon 
751       goto (10,20,30,30,40) ipot
752 C----------------------- LJ potential ---------------------------------
753    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
754      &   (sigma0(i),i=1,ntyp)
755       if (lprint) then
756         write (iout,'(/a/)') 'Parameters of the LJ potential:'
757         write (iout,'(a/)') 'The epsilon array:'
758         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
759         write (iout,'(/a)') 'One-body parameters:'
760         write (iout,'(a,4x,a)') 'residue','sigma'
761         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
762       endif
763       goto 50
764 C----------------------- LJK potential --------------------------------
765    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
766      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
767       if (lprint) then
768         write (iout,'(/a/)') 'Parameters of the LJK potential:'
769         write (iout,'(a/)') 'The epsilon array:'
770         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
771         write (iout,'(/a)') 'One-body parameters:'
772         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
773         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
774      &        i=1,ntyp)
775       endif
776       goto 50
777 C---------------------- GB or BP potential -----------------------------
778    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
779      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
780      &  (alp(i),i=1,ntyp)
781 C For the GB potential convert sigma'**2 into chi'
782       if (ipot.eq.4) then
783         do i=1,ntyp
784           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
785         enddo
786       endif
787       if (lprint) then
788         write (iout,'(/a/)') 'Parameters of the BP potential:'
789         write (iout,'(a/)') 'The epsilon array:'
790         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
791         write (iout,'(/a)') 'One-body parameters:'
792         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
793      &       '    chip  ','    alph  '
794         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
795      &                     chip(i),alp(i),i=1,ntyp)
796       endif
797       goto 50
798 C--------------------- GBV potential -----------------------------------
799    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
800      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
801      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
802       if (lprint) then
803         write (iout,'(/a/)') 'Parameters of the GBV potential:'
804         write (iout,'(a/)') 'The epsilon array:'
805         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
806         write (iout,'(/a)') 'One-body parameters:'
807         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
808      &      's||/s_|_^2','    chip  ','    alph  '
809         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
810      &           sigii(i),chip(i),alp(i),i=1,ntyp)
811       endif
812    50 continue
813       close (isidep)
814 C-----------------------------------------------------------------------
815 C Calculate the "working" parameters of SC interactions.
816       do i=2,ntyp
817         do j=1,i-1
818           eps(i,j)=eps(j,i)
819         enddo
820       enddo
821       do i=1,ntyp
822         do j=i,ntyp
823           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
824           sigma(j,i)=sigma(i,j)
825           rs0(i,j)=dwa16*sigma(i,j)
826           rs0(j,i)=rs0(i,j)
827         enddo
828       enddo
829       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
830      & 'Working parameters of the SC interactions:',
831      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
832      & '  chi1   ','   chi2   ' 
833       do i=1,ntyp
834         do j=i,ntyp
835           epsij=eps(i,j)
836           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
837             rrij=sigma(i,j)
838           else
839             rrij=rr0(i)+rr0(j)
840           endif
841           r0(i,j)=rrij
842           r0(j,i)=rrij
843           rrij=rrij**expon
844           epsij=eps(i,j)
845           sigeps=dsign(1.0D0,epsij)
846           epsij=dabs(epsij)
847           aa(i,j)=epsij*rrij*rrij
848           bb(i,j)=-sigeps*epsij*rrij
849           aa(j,i)=aa(i,j)
850           bb(j,i)=bb(i,j)
851           if (ipot.gt.2) then
852             sigt1sq=sigma0(i)**2
853             sigt2sq=sigma0(j)**2
854             sigii1=sigii(i)
855             sigii2=sigii(j)
856             ratsig1=sigt2sq/sigt1sq
857             ratsig2=1.0D0/ratsig1
858             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
859             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
860             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
861           else
862             rsum_max=sigma(i,j)
863           endif
864 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
865             sigmaii(i,j)=rsum_max
866             sigmaii(j,i)=rsum_max 
867 c         else
868 c           sigmaii(i,j)=r0(i,j)
869 c           sigmaii(j,i)=r0(i,j)
870 c         endif
871 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
872           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
873             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
874             augm(i,j)=epsij*r_augm**(2*expon)
875 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
876             augm(j,i)=augm(i,j)
877           else
878             augm(i,j)=0.0D0
879             augm(j,i)=0.0D0
880           endif
881           if (lprint) then
882             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
883      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
884      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
885           endif
886         enddo
887       enddo
888 #ifdef OLDSCP
889 C
890 C Define the SC-p interaction constants (hard-coded; old style)
891 C
892       do i=1,20
893 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
894 C helix formation)
895 c       aad(i,1)=0.3D0*4.0D0**12
896 C Following line for constants currently implemented
897 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
898         aad(i,1)=1.5D0*4.0D0**12
899 c       aad(i,1)=0.17D0*5.6D0**12
900         aad(i,2)=aad(i,1)
901 C "Soft" SC-p repulsion
902         bad(i,1)=0.0D0
903 C Following line for constants currently implemented
904 c       aad(i,1)=0.3D0*4.0D0**6
905 C "Hard" SC-p repulsion
906         bad(i,1)=3.0D0*4.0D0**6
907 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
908         bad(i,2)=bad(i,1)
909 c       aad(i,1)=0.0D0
910 c       aad(i,2)=0.0D0
911 c       bad(i,1)=1228.8D0
912 c       bad(i,2)=1228.8D0
913       enddo
914 #else
915 C
916 C 8/9/01 Read the SC-p interaction constants from file
917 C
918       do i=1,ntyp
919         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
920       enddo
921       do i=1,ntyp
922         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
923         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
924         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
925         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
926       enddo
927
928       if (lprint) then
929         write (iout,*) "Parameters of SC-p interactions:"
930         do i=1,20
931           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
932      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
933         enddo
934       endif
935 #endif
936 C
937 C Define the constants of the disulfide bridge
938 C
939       ebr=-5.50D0
940 c
941 c Old arbitrary potential - commented out.
942 c
943 c      dbr= 4.20D0
944 c      fbr= 3.30D0
945 c
946 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
947 c energy surface of diethyl disulfide.
948 c A. Liwo and U. Kozlowska, 11/24/03
949 c
950       D0CM = 3.78d0
951       AKCM = 15.1d0
952       AKTH = 11.0d0
953       AKCT = 12.0d0
954       V1SS =-1.08d0
955       V2SS = 7.61d0
956       V3SS = 13.7d0
957 c      akcm=0.0d0
958 c      akth=0.0d0
959 c      akct=0.0d0
960 c      v1ss=0.0d0
961 c      v2ss=0.0d0
962 c      v3ss=0.0d0
963       
964       if(me.eq.king) then
965       write (iout,'(/a)') "Disulfide bridge parameters:"
966       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
967       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
968       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
969       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
970      &  ' v3ss:',v3ss
971       endif
972       return
973   111 write (iout,*) "Error reading bending energy parameters."
974       goto 999
975   112 write (iout,*) "Error reading rotamer energy parameters."
976       goto 999
977   113 write (iout,*) "Error reading torsional energy parameters."
978       goto 999
979   114 write (iout,*) "Error reading double torsional energy parameters."
980       goto 999
981   115 write (iout,*) 
982      &  "Error reading cumulant (multibody energy) parameters."
983       goto 999
984   116 write (iout,*) "Error reading electrostatic energy parameters."
985       goto 999
986   117 write (iout,*) "Error reading side chain interaction parameters."
987       goto 999
988   118 write (iout,*) "Error reading SCp interaction parameters."
989       goto 999
990   119 write (iout,*) "Error reading SCCOR parameters"
991   999 continue
992 #ifdef MPI
993       call MPI_Finalize(Ierror)
994 #endif
995       stop
996       return
997       end
998
999
1000       subroutine getenv_loc(var, val)
1001       character(*) var, val
1002
1003 #ifdef WINIFL
1004       character(2000) line
1005       external ilen
1006
1007       open (196,file='env',status='old',readonly,shared)
1008       iread=0
1009 c      write(*,*)'looking for ',var
1010 10    read(196,*,err=11,end=11)line
1011       iread=index(line,var)
1012 c      write(*,*)iread,' ',var,' ',line
1013       if (iread.eq.0) go to 10 
1014 c      write(*,*)'---> ',line
1015 11    continue
1016       if(iread.eq.0) then
1017 c       write(*,*)'CHUJ'
1018        val=''
1019       else
1020        iread=iread+ilen(var)+1
1021        read (line(iread:),*,err=12,end=12) val
1022 c       write(*,*)'OK: ',var,' = ',val
1023       endif
1024       close(196)
1025       return
1026 12    val=''
1027       close(196)
1028 #elif (defined CRAY)
1029       integer lennam,lenval,ierror
1030 c
1031 c        getenv using a POSIX call, useful on the T3D
1032 c        Sept 1996, comment out error check on advice of H. Pritchard
1033 c
1034       lennam = len(var)
1035       if(lennam.le.0) stop '--error calling getenv--'
1036       call pxfgetenv(var,lennam,val,lenval,ierror)
1037 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1038 #else
1039       call getenv(var,val)
1040 #endif
1041
1042       return
1043       end