574552cdf0ff69012b920b1c3ee8bfd0da988e9a
[unres.git] / 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       include 'COMMON.CONTROL'
30       include 'COMMON.SHIELD'
31       character*1 t1,t2,t3
32       character*1 onelett(4) /"G","A","P","D"/
33       character*1 toronelet(-2:2) /"p","a","G","A","P"/
34       logical lprint,LaTeX
35       dimension blower(3,3,maxlob)
36 C      dimension b(13)
37       character*3 lancuch,ucase
38 C
39 C For printing parameters after they are read set the following in the UNRES
40 C C-shell script:
41 C
42 C setenv PRINT_PARM YES
43 C
44 C To print parameters in LaTeX format rather than as ASCII tables:
45 C
46 C setenv LATEX YES
47 C
48       call getenv_loc("PRINT_PARM",lancuch)
49       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
50       call getenv_loc("LATEX",lancuch)
51       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
52 C
53       dwa16=2.0d0**(1.0d0/6.0d0)
54       itypro=20
55 C Assign virtual-bond length
56       vbl=3.8D0
57       vblinv=1.0D0/vbl
58       vblinv2=vblinv*vblinv
59 c
60 c Read the virtual-bond parameters, masses, and moments of inertia
61 c and Stokes' radii of the peptide group and side chains
62 c
63 #ifdef CRYST_BOND
64       read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
65       do i=1,ntyp
66         nbondterm(i)=1
67         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
68         dsc(i) = vbldsc0(1,i)
69         if (i.eq.10) then
70           dsc_inv(i)=0.0D0
71         else
72           dsc_inv(i)=1.0D0/dsc(i)
73         endif
74       enddo
75 #else
76       read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
77       do i=1,ntyp
78       print *,i
79         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
80      &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
81         dsc(i) = vbldsc0(1,i)
82         if (i.eq.10) then
83           dsc_inv(i)=0.0D0
84         else
85           dsc_inv(i)=1.0D0/dsc(i)
86         endif
87       enddo
88 #endif
89       if (lprint) then
90         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
91         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
92      &   'inertia','Pstok'
93         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
94         do i=1,ntyp
95           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
96      &      vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
97           do j=2,nbondterm(i)
98             write (iout,'(13x,3f10.5)')
99      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
100           enddo
101         enddo
102       endif
103 C reading lipid parameters
104        read(iliptranpar,*) pepliptran
105        do i=1,ntyp
106        read(iliptranpar,*) liptranene(i)
107        enddo
108        close(iliptranpar)
109 #ifdef CRYST_THETA
110 C
111 C Read the parameters of the probability distribution/energy expression 
112 C of the virtual-bond valence angles theta
113 C
114       do i=1,ntyp
115         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
116      &    (bthet(j,i,1,1),j=1,2)
117         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
118         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
119         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
120         sigc0(i)=sigc0(i)**2
121       enddo
122       do i=1,ntyp
123       athet(1,i,1,-1)=athet(1,i,1,1)
124       athet(2,i,1,-1)=athet(2,i,1,1)
125       bthet(1,i,1,-1)=-bthet(1,i,1,1)
126       bthet(2,i,1,-1)=-bthet(2,i,1,1)
127       athet(1,i,-1,1)=-athet(1,i,1,1)
128       athet(2,i,-1,1)=-athet(2,i,1,1)
129       bthet(1,i,-1,1)=bthet(1,i,1,1)
130       bthet(2,i,-1,1)=bthet(2,i,1,1)
131       enddo
132       do i=-ntyp,-1
133       a0thet(i)=a0thet(-i)
134       athet(1,i,-1,-1)=athet(1,-i,1,1)
135       athet(2,i,-1,-1)=-athet(2,-i,1,1)
136       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
137       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
138       athet(1,i,-1,1)=athet(1,-i,1,1)
139       athet(2,i,-1,1)=-athet(2,-i,1,1)
140       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
141       bthet(2,i,-1,1)=bthet(2,-i,1,1)
142       athet(1,i,1,-1)=-athet(1,-i,1,1)
143       athet(2,i,1,-1)=athet(2,-i,1,1)
144       bthet(1,i,1,-1)=bthet(1,-i,1,1)
145       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
146       theta0(i)=theta0(-i)
147       sig0(i)=sig0(-i)
148       sigc0(i)=sigc0(-i)
149        do j=0,3
150         polthet(j,i)=polthet(j,-i)
151        enddo
152        do j=1,3
153          gthet(j,i)=gthet(j,-i)
154        enddo
155       enddo
156
157       close (ithep)
158       if (lprint) then
159       if (.not.LaTeX) then
160         write (iout,'(a)') 
161      &    'Parameters of the virtual-bond valence angles:'
162         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
163      & '    ATHETA0   ','         A1   ','        A2    ',
164      & '        B1    ','         B2   '        
165         do i=1,ntyp
166           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
167      &        a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
168         enddo
169         write (iout,'(/a/9x,5a/79(1h-))') 
170      & 'Parameters of the expression for sigma(theta_c):',
171      & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
172      & '     ALPH3    ','    SIGMA0C   '        
173         do i=1,ntyp
174           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
175      &      (polthet(j,i),j=0,3),sigc0(i) 
176         enddo
177         write (iout,'(/a/9x,5a/79(1h-))') 
178      & 'Parameters of the second gaussian:',
179      & '    THETA0    ','     SIGMA0   ','        G1    ',
180      & '        G2    ','         G3   '        
181         do i=1,ntyp
182           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
183      &       sig0(i),(gthet(j,i),j=1,3)
184         enddo
185        else
186         write (iout,'(a)') 
187      &    'Parameters of the virtual-bond valence angles:'
188         write (iout,'(/a/9x,5a/79(1h-))') 
189      & 'Coefficients of expansion',
190      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
191      & '   b1*10^1    ','    b2*10^1   '        
192         do i=1,ntyp
193           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
194      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
195      &        (10*bthet(j,i,1,1),j=1,2)
196         enddo
197         write (iout,'(/a/9x,5a/79(1h-))') 
198      & 'Parameters of the expression for sigma(theta_c):',
199      & ' alpha0       ','  alph1       ',' alph2        ',
200      & ' alhp3        ','   sigma0c    '        
201         do i=1,ntyp
202           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
203      &      (polthet(j,i),j=0,3),sigc0(i) 
204         enddo
205         write (iout,'(/a/9x,5a/79(1h-))') 
206      & 'Parameters of the second gaussian:',
207      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
208      & '        G2    ','   G3*10^1    '        
209         do i=1,ntyp
210           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
211      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
212         enddo
213       endif
214       endif
215 #else 
216
217 C Read the parameters of Utheta determined from ab initio surfaces
218 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
219 C
220       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
221      &  ntheterm3,nsingle,ndouble
222       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
223       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
224       do i=-ntyp1,-1
225         ithetyp(i)=-ithetyp(-i)
226       enddo
227       do iblock=1,2
228       do i=-maxthetyp,maxthetyp
229         do j=-maxthetyp,maxthetyp
230           do k=-maxthetyp,maxthetyp
231             aa0thet(i,j,k,iblock)=0.0d0
232             do l=1,ntheterm
233               aathet(l,i,j,k,iblock)=0.0d0
234             enddo
235             do l=1,ntheterm2
236               do m=1,nsingle
237                 bbthet(m,l,i,j,k,iblock)=0.0d0
238                 ccthet(m,l,i,j,k,iblock)=0.0d0
239                 ddthet(m,l,i,j,k,iblock)=0.0d0
240                 eethet(m,l,i,j,k,iblock)=0.0d0
241               enddo
242             enddo
243             do l=1,ntheterm3
244               do m=1,ndouble
245                 do mm=1,ndouble
246                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
247                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
248                 enddo
249               enddo
250             enddo
251           enddo
252         enddo
253       enddo
254       enddo
255 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
256       do iblock=1,2 
257 c VAR:ntethtyp is type of theta potentials type currently 0=glycine 
258 c VAR:1=non-glicyne non-proline 2=proline
259 c VAR:negative values for D-aminoacid
260       do i=0,nthetyp
261         do j=-nthetyp,nthetyp
262           do k=-nthetyp,nthetyp
263             read (ithep,'(6a)',end=111,err=111) res1
264             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
265 c VAR: aa0thet is variable describing the average value of Foureir
266 c VAR: expansion series
267 c VAR: aathet is foureir expansion in theta/2 angle for full formula
268 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
269 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
270             read (ithep,*,end=111,err=111)
271      &(aathet(l,i,j,k,iblock),l=1,ntheterm)
272             read (ithep,*,end=111,err=111)
273      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
274      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
275      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
276      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
277      &        ll=1,ntheterm2)
278             read (ithep,*,end=111,err=111)
279      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
280      &         ffthet(lll,llll,ll,i,j,k,iblock),
281      &         ggthet(llll,lll,ll,i,j,k,iblock),
282      &         ggthet(lll,llll,ll,i,j,k,iblock),
283      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
284           enddo
285         enddo
286       enddo
287 C
288 C For dummy ends assign glycine-type coefficients of theta-only terms; the
289 C coefficients of theta-and-gamma-dependent terms are zero.
290 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
291 C RECOMENTDED AFTER VERSION 3.3)
292 c      do i=1,nthetyp
293 c        do j=1,nthetyp
294 c          do l=1,ntheterm
295 c            aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
296 c            aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
297 c          enddo
298 c          aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
299 c          aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
300 c        enddo
301 c        do l=1,ntheterm
302 c          aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
303 c        enddo
304 c        aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
305 c      enddo
306 c      enddo
307 C AND COMMENT THE LOOPS BELOW
308       do i=1,nthetyp
309         do j=1,nthetyp
310           do l=1,ntheterm
311             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
312             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
313           enddo
314           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
315           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
316         enddo
317         do l=1,ntheterm
318           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
319         enddo
320         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
321       enddo
322       enddo
323 C TILL HERE
324 C Substitution for D aminoacids from symmetry.
325       do iblock=1,2
326       do i=-nthetyp,0
327         do j=-nthetyp,nthetyp
328           do k=-nthetyp,nthetyp
329            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
330            do l=1,ntheterm
331            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock) 
332            enddo
333            do ll=1,ntheterm2
334             do lll=1,nsingle
335             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
336             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
337             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
338             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
339             enddo
340           enddo
341           do ll=1,ntheterm3
342            do lll=2,ndouble
343             do llll=1,lll-1
344             ffthet(llll,lll,ll,i,j,k,iblock)=
345      &      ffthet(llll,lll,ll,-i,-j,-k,iblock) 
346             ffthet(lll,llll,ll,i,j,k,iblock)=
347      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
348             ggthet(llll,lll,ll,i,j,k,iblock)=
349      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
350             ggthet(lll,llll,ll,i,j,k,iblock)=
351      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)      
352             enddo !ll
353            enddo  !lll  
354           enddo   !llll
355          enddo    !k
356         enddo     !j
357        enddo      !i
358       enddo       !iblock
359 C
360 C Control printout of the coefficients of virtual-bond-angle potentials
361 C
362       if (lprint) then
363         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
364         do i=1,nthetyp+1
365           do j=1,nthetyp+1
366             do k=1,nthetyp+1
367               write (iout,'(//4a)') 
368      &         'Type ',onelett(i),onelett(j),onelett(k) 
369               write (iout,'(//a,10x,a)') " l","a[l]"
370               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
371               write (iout,'(i2,1pe15.5)')
372      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
373             do l=1,ntheterm2
374               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') 
375      &          "b",l,"c",l,"d",l,"e",l
376               do m=1,nsingle
377                 write (iout,'(i2,4(1pe15.5))') m,
378      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
379      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
380               enddo
381             enddo
382             do l=1,ntheterm3
383               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
384      &          "f+",l,"f-",l,"g+",l,"g-",l
385               do m=2,ndouble
386                 do n=1,m-1
387                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
388      &              ffthet(n,m,l,i,j,k,iblock),
389      &              ffthet(m,n,l,i,j,k,iblock),
390      &              ggthet(n,m,l,i,j,k,iblock),
391      &              ggthet(m,n,l,i,j,k,iblock)
392                 enddo
393               enddo
394             enddo
395           enddo
396         enddo
397       enddo
398       call flush(iout)
399       endif
400       write (2,*) "Start reading THETA_PDB",ithep_pdb
401       do i=1,ntyp
402 c      write (2,*) 'i=',i
403         read (ithep_pdb,*,err=111,end=111)
404      &     a0thet(i),(athet(j,i,1,1),j=1,2),
405      &    (bthet(j,i,1,1),j=1,2)
406         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
407         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
408         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
409         sigc0(i)=sigc0(i)**2
410       enddo
411       do i=1,ntyp
412       athet(1,i,1,-1)=athet(1,i,1,1)
413       athet(2,i,1,-1)=athet(2,i,1,1)
414       bthet(1,i,1,-1)=-bthet(1,i,1,1)
415       bthet(2,i,1,-1)=-bthet(2,i,1,1)
416       athet(1,i,-1,1)=-athet(1,i,1,1)
417       athet(2,i,-1,1)=-athet(2,i,1,1)
418       bthet(1,i,-1,1)=bthet(1,i,1,1)
419       bthet(2,i,-1,1)=bthet(2,i,1,1)
420       enddo
421       do i=-ntyp,-1
422       a0thet(i)=a0thet(-i)
423       athet(1,i,-1,-1)=athet(1,-i,1,1)
424       athet(2,i,-1,-1)=-athet(2,-i,1,1)
425       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
426       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
427       athet(1,i,-1,1)=athet(1,-i,1,1)
428       athet(2,i,-1,1)=-athet(2,-i,1,1)
429       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
430       bthet(2,i,-1,1)=bthet(2,-i,1,1)
431       athet(1,i,1,-1)=-athet(1,-i,1,1)
432       athet(2,i,1,-1)=athet(2,-i,1,1)
433       bthet(1,i,1,-1)=bthet(1,-i,1,1)
434       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
435       theta0(i)=theta0(-i)
436       sig0(i)=sig0(-i)
437       sigc0(i)=sigc0(-i)
438        do j=0,3
439         polthet(j,i)=polthet(j,-i)
440        enddo
441        do j=1,3
442          gthet(j,i)=gthet(j,-i)
443        enddo
444       enddo
445       write (2,*) "End reading THETA_PDB"
446       close (ithep_pdb)
447 #endif
448       close(ithep)
449 #ifdef CRYST_SC
450 C
451 C Read the parameters of the probability distribution/energy expression
452 C of the side chains.
453 C
454       do i=1,ntyp
455         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
456         if (i.eq.10) then
457           dsc_inv(i)=0.0D0
458         else
459           dsc_inv(i)=1.0D0/dsc(i)
460         endif
461         if (i.ne.10) then
462         do j=1,nlob(i)
463           do k=1,3
464             do l=1,3
465               blower(l,k,j)=0.0D0
466             enddo
467           enddo
468         enddo  
469         bsc(1,i)=0.0D0
470         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
471      &    ((blower(k,l,1),l=1,k),k=1,3)
472         censc(1,1,-i)=censc(1,1,i)
473         censc(2,1,-i)=censc(2,1,i)
474         censc(3,1,-i)=-censc(3,1,i)
475         do j=2,nlob(i)
476           read (irotam,*,end=112,err=112) bsc(j,i)
477           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
478      &                                 ((blower(k,l,j),l=1,k),k=1,3)
479         censc(1,j,-i)=censc(1,j,i)
480         censc(2,j,-i)=censc(2,j,i)
481         censc(3,j,-i)=-censc(3,j,i)
482 C BSC is amplitude of Gaussian
483         enddo
484         do j=1,nlob(i)
485           do k=1,3
486             do l=1,k
487               akl=0.0D0
488               do m=1,3
489                 akl=akl+blower(k,m,j)*blower(l,m,j)
490               enddo
491               gaussc(k,l,j,i)=akl
492               gaussc(l,k,j,i)=akl
493              if (((k.eq.3).and.(l.ne.3))
494      &        .or.((l.eq.3).and.(k.ne.3))) then
495                 gaussc(k,l,j,-i)=-akl
496                 gaussc(l,k,j,-i)=-akl
497               else
498                 gaussc(k,l,j,-i)=akl
499                 gaussc(l,k,j,-i)=akl
500               endif
501             enddo
502           enddo 
503         enddo
504         endif
505       enddo
506       close (irotam)
507       if (lprint) then
508         write (iout,'(/a)') 'Parameters of side-chain local geometry'
509         do i=1,ntyp
510           nlobi=nlob(i)
511           if (nlobi.gt.0) then
512             if (LaTeX) then
513               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
514      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
515                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
516      &                             'log h',(bsc(j,i),j=1,nlobi)
517                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
518      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
519               do k=1,3
520                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
521      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
522               enddo
523             else
524               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
525               write (iout,'(a,f10.4,4(16x,f10.4))')
526      &                             'Center  ',(bsc(j,i),j=1,nlobi)
527               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
528      &           j=1,nlobi)
529               write (iout,'(a)')
530             endif
531           endif
532         enddo
533       endif
534 #else
535
536 C Read scrot parameters for potentials determined from all-atom AM1 calculations
537 C added by Urszula Kozlowska 07/11/2007
538 C
539       do i=1,ntyp
540         read (irotam,*,end=112,err=112) 
541        if (i.eq.10) then 
542          read (irotam,*,end=112,err=112) 
543        else
544          do j=1,65
545            read(irotam,*,end=112,err=112) sc_parmin(j,i)
546          enddo  
547        endif
548       enddo
549 C
550 C Read the parameters of the probability distribution/energy expression
551 C of the side chains.
552 C
553       write (2,*) "Start reading ROTAM_PDB"
554       do i=1,ntyp
555         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
556         if (i.eq.10) then
557           dsc_inv(i)=0.0D0
558         else
559           dsc_inv(i)=1.0D0/dsc(i)
560         endif
561         if (i.ne.10) then
562         do j=1,nlob(i)
563           do k=1,3
564             do l=1,3
565               blower(l,k,j)=0.0D0
566             enddo
567           enddo
568         enddo
569         bsc(1,i)=0.0D0
570         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
571      &    ((blower(k,l,1),l=1,k),k=1,3)
572         do j=2,nlob(i)
573           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
574           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
575      &                                 ((blower(k,l,j),l=1,k),k=1,3)
576         enddo
577         do j=1,nlob(i)
578           do k=1,3
579             do l=1,k
580               akl=0.0D0
581               do m=1,3
582                 akl=akl+blower(k,m,j)*blower(l,m,j)
583               enddo
584               gaussc(k,l,j,i)=akl
585               gaussc(l,k,j,i)=akl
586             enddo
587           enddo
588         enddo
589         endif
590       enddo
591       close (irotam_pdb)
592       write (2,*) "End reading ROTAM_PDB"
593 #endif
594       close(irotam)
595
596 #ifdef CRYST_TOR
597 C
598 C Read torsional parameters in old format
599 C
600       read (itorp,*,end=113,err=113) ntortyp,nterm_old
601       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
602       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
603       do i=1,ntortyp
604         do j=1,ntortyp
605           read (itorp,'(a)')
606           do k=1,nterm_old
607             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i) 
608           enddo
609         enddo
610       enddo
611       close (itorp)
612       if (lprint) then
613         write (iout,'(/a/)') 'Torsional constants:'
614         do i=1,ntortyp
615           do j=1,ntortyp
616             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
617             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
618           enddo
619         enddo
620       endif
621 #else
622 C
623 C Read torsional parameters
624 C
625       read (itorp,*,end=113,err=113) ntortyp
626       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
627       do iblock=1,2
628       do i=-ntyp,-1
629        itortyp(i)=-itortyp(-i)
630       enddo
631       write (iout,*) 'ntortyp',ntortyp
632       do i=0,ntortyp-1
633         do j=-ntortyp+1,ntortyp-1
634           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
635      &          nlor(i,j,iblock)
636           nterm(-i,-j,iblock)=nterm(i,j,iblock)
637           nlor(-i,-j,iblock)=nlor(i,j,iblock)
638           v0ij=0.0d0
639           si=-1.0d0
640           do k=1,nterm(i,j,iblock)
641             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
642      &      v2(k,i,j,iblock)
643             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
644             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
645             v0ij=v0ij+si*v1(k,i,j,iblock)
646             si=-si
647 c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
648 c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
649 c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
650           enddo
651           do k=1,nlor(i,j,iblock)
652             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
653      &        vlor2(k,i,j),vlor3(k,i,j)
654             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
655           enddo
656           v0(i,j,iblock)=v0ij
657           v0(-i,-j,iblock)=v0ij
658         enddo
659       enddo
660       enddo
661       close (itorp)
662       if (lprint) then
663         write (iout,'(/a/)') 'Torsional constants:'
664         do i=1,ntortyp
665           do j=1,ntortyp
666             write (iout,*) 'ityp',i,' jtyp',j
667             write (iout,*) 'Fourier constants'
668             do k=1,nterm(i,j,iblock)
669               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
670      &        v2(k,i,j,iblock)
671             enddo
672             write (iout,*) 'Lorenz constants'
673             do k=1,nlor(i,j,iblock)
674               write (iout,'(3(1pe15.5))')
675      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
676             enddo
677           enddo
678         enddo
679       endif
680
681 C
682 C 6/23/01 Read parameters for double torsionals
683 C
684       do iblock=1,2
685       do i=0,ntortyp-1
686         do j=-ntortyp+1,ntortyp-1
687           do k=-ntortyp+1,ntortyp-1
688             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
689 c              write (iout,*) "OK onelett",
690 c     &         i,j,k,t1,t2,t3
691
692             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
693      &        .or. t3.ne.toronelet(k)) then
694               write (iout,*) "Error in double torsional parameter file",
695      &         i,j,k,t1,t2,t3
696 #ifdef MPI
697               call MPI_Finalize(Ierror)
698 #endif
699                stop "Error in double torsional parameter file"
700             endif
701            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
702      &         ntermd_2(i,j,k,iblock)
703             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
704             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
705             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
706      &         ntermd_1(i,j,k,iblock))
707             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
708      &         ntermd_1(i,j,k,iblock))
709             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
710      &         ntermd_1(i,j,k,iblock))
711             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
712      &         ntermd_1(i,j,k,iblock))
713 C Martix of D parameters for one dimesional foureir series
714             do l=1,ntermd_1(i,j,k,iblock)
715              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
716              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
717              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
718              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
719 c            write(iout,*) "whcodze" ,
720 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
721             enddo
722             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
723      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
724      &         v2s(m,l,i,j,k,iblock),
725      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
726 C Martix of D parameters for two dimesional fourier series
727             do l=1,ntermd_2(i,j,k,iblock)
728              do m=1,l-1
729              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
730              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
731              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
732              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
733              enddo!m
734             enddo!l
735           enddo!k
736         enddo!j
737       enddo!i
738       enddo!iblock
739       if (lprint) then
740       write (iout,*)
741       write (iout,*) 'Constants for double torsionals'
742       do iblock=1,2
743       do i=0,ntortyp-1
744         do j=-ntortyp+1,ntortyp-1
745           do k=-ntortyp+1,ntortyp-1
746             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
747      &        ' nsingle',ntermd_1(i,j,k,iblock),
748      &        ' ndouble',ntermd_2(i,j,k,iblock)
749             write (iout,*)
750             write (iout,*) 'Single angles:'
751             do l=1,ntermd_1(i,j,k,iblock)
752               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
753      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
754      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
755      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
756             enddo
757             write (iout,*)
758             write (iout,*) 'Pairs of angles:'
759             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
760             do l=1,ntermd_2(i,j,k,iblock)
761               write (iout,'(i5,20f10.5)')
762      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
763             enddo
764             write (iout,*)
765             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
766             do l=1,ntermd_2(i,j,k,iblock)
767               write (iout,'(i5,20f10.5)')
768      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
769      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
770             enddo
771             write (iout,*)
772           enddo
773         enddo
774       enddo
775       enddo
776       endif
777 #endif
778 C Read of Side-chain backbone correlation parameters
779 C Modified 11 May 2012 by Adasko
780 CCC
781 C
782       read (isccor,*,end=119,err=119) nsccortyp
783 #ifdef SCCORPDB
784       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
785       do i=-ntyp,-1
786         isccortyp(i)=-isccortyp(-i)
787       enddo
788       iscprol=isccortyp(20)
789 c      write (iout,*) 'ntortyp',ntortyp
790       maxinter=3
791 cc maxinter is maximum interaction sites
792       do l=1,maxinter
793       do i=1,nsccortyp
794         do j=1,nsccortyp
795           read (isccor,*,end=119,err=119)
796      &nterm_sccor(i,j),nlor_sccor(i,j)
797           v0ijsccor=0.0d0
798           v0ijsccor1=0.0d0
799           v0ijsccor2=0.0d0
800           v0ijsccor3=0.0d0
801           si=-1.0d0
802           nterm_sccor(-i,j)=nterm_sccor(i,j)
803           nterm_sccor(-i,-j)=nterm_sccor(i,j)
804           nterm_sccor(i,-j)=nterm_sccor(i,j)
805           do k=1,nterm_sccor(i,j)
806             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
807      &    ,v2sccor(k,l,i,j)
808             if (j.eq.iscprol) then
809              if (i.eq.isccortyp(10)) then
810              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
811              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
812              else
813              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
814      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
815              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
816      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
817              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
818              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
819              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
820              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
821              endif
822             else
823              if (i.eq.isccortyp(10)) then
824              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
825              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
826              else
827                if (j.eq.isccortyp(10)) then
828              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
829              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
830                else
831              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
832              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
833              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
834              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
835              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
836              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
837                 endif
838                endif
839             endif
840             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
841             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
842             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
843             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
844             si=-si
845           enddo
846           do k=1,nlor_sccor(i,j)
847             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
848      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
849             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
850      &(1+vlor3sccor(k,i,j)**2)
851           enddo
852           v0sccor(l,i,j)=v0ijsccor
853           v0sccor(l,-i,j)=v0ijsccor1
854           v0sccor(l,i,-j)=v0ijsccor2
855           v0sccor(l,-i,-j)=v0ijsccor3         
856         enddo
857       enddo
858       enddo
859       close (isccor)
860 #else
861       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
862 c      write (iout,*) 'ntortyp',ntortyp
863       maxinter=3
864 cc maxinter is maximum interaction sites
865       do l=1,maxinter
866       do i=1,nsccortyp
867         do j=1,nsccortyp
868           read (isccor,*,end=119,err=119)
869      & nterm_sccor(i,j),nlor_sccor(i,j)
870           v0ijsccor=0.0d0
871           si=-1.0d0
872
873           do k=1,nterm_sccor(i,j)
874             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
875      &    ,v2sccor(k,l,i,j)
876             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
877             si=-si
878           enddo
879           do k=1,nlor_sccor(i,j)
880             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
881      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
882             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
883      &(1+vlor3sccor(k,i,j)**2)
884           enddo
885           v0sccor(i,j,iblock)=v0ijsccor
886         enddo
887       enddo
888       enddo
889       close (isccor)
890
891 #endif      
892       if (lprint) then
893         write (iout,'(/a/)') 'Torsional constants:'
894         do i=1,nsccortyp
895           do j=1,nsccortyp
896             write (iout,*) 'ityp',i,' jtyp',j
897             write (iout,*) 'Fourier constants'
898             do k=1,nterm_sccor(i,j)
899       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
900             enddo
901             write (iout,*) 'Lorenz constants'
902             do k=1,nlor_sccor(i,j)
903               write (iout,'(3(1pe15.5))')
904      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
905             enddo
906           enddo
907         enddo
908       endif
909
910 C
911 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
912 C         interaction energy of the Gly, Ala, and Pro prototypes.
913 C
914       if (lprint) then
915         write (iout,*)
916         write (iout,*) "Coefficients of the cumulants"
917       endif
918       read (ifourier,*) nloctyp
919
920       do i=0,nloctyp-1
921         read (ifourier,*,end=115,err=115)
922         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
923 #ifdef NEWCORR
924         read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
925         read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
926         read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
927         read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
928         read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
929 #endif 
930         if (lprint) then
931         write (iout,*) 'Type',i
932         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
933         endif
934 c        B1(1,i)  = b(3)
935 c        B1(2,i)  = b(5)
936 c        B1(1,-i) = b(3)
937 c        B1(2,-i) = -b(5)
938 c        b1(1,i)=0.0d0
939 c        b1(2,i)=0.0d0
940 c        B1tilde(1,i) = b(3)
941 c        B1tilde(2,i) =-b(5)
942 c        B1tilde(1,-i) =-b(3)
943 c        B1tilde(2,-i) =b(5)
944 c        b1tilde(1,i)=0.0d0
945 c        b1tilde(2,i)=0.0d0
946 c        B2(1,i)  = b(2)
947 c        B2(2,i)  = b(4)
948 c        B2(1,-i)  =b(2)
949 c        B2(2,-i)  =-b(4)
950         B1tilde(1,i) = b(3,i)
951         B1tilde(2,i) =-b(5,i)
952         B1tilde(1,-i) =-b(3,i)
953         B1tilde(2,-i) =b(5,i)
954         b1tilde(1,i)=0.0d0
955         b1tilde(2,i)=0.0d0
956         B2(1,i)  = b(2,i)
957         B2(2,i)  = b(4,i)
958         B2(1,-i)  =b(2,i)
959         B2(2,-i)  =-b(4,i)
960
961 c        b2(1,i)=0.0d0
962 c        b2(2,i)=0.0d0
963         CC(1,1,i)= b(7,i)
964         CC(2,2,i)=-b(7,i)
965         CC(2,1,i)= b(9,i)
966         CC(1,2,i)= b(9,i)
967         CC(1,1,-i)= b(7,i)
968         CC(2,2,-i)=-b(7,i)
969         CC(2,1,-i)=-b(9,i)
970         CC(1,2,-i)=-b(9,i)
971 c        CC(1,1,i)=0.0d0
972 c        CC(2,2,i)=0.0d0
973 c        CC(2,1,i)=0.0d0
974 c        CC(1,2,i)=0.0d0
975         Ctilde(1,1,i)=b(7,i)
976         Ctilde(1,2,i)=b(9,i)
977         Ctilde(2,1,i)=-b(9,i)
978         Ctilde(2,2,i)=b(7,i)
979         Ctilde(1,1,-i)=b(7,i)
980         Ctilde(1,2,-i)=-b(9,i)
981         Ctilde(2,1,-i)=b(9,i)
982         Ctilde(2,2,-i)=b(7,i)
983
984 c        Ctilde(1,1,i)=0.0d0
985 c        Ctilde(1,2,i)=0.0d0
986 c        Ctilde(2,1,i)=0.0d0
987 c        Ctilde(2,2,i)=0.0d0
988         DD(1,1,i)= b(6,i)
989         DD(2,2,i)=-b(6,i)
990         DD(2,1,i)= b(8,i)
991         DD(1,2,i)= b(8,i)
992         DD(1,1,-i)= b(6,i)
993         DD(2,2,-i)=-b(6,i)
994         DD(2,1,-i)=-b(8,i)
995         DD(1,2,-i)=-b(8,i)
996 c        DD(1,1,i)=0.0d0
997 c        DD(2,2,i)=0.0d0
998 c        DD(2,1,i)=0.0d0
999 c        DD(1,2,i)=0.0d0
1000         Dtilde(1,1,i)=b(6,i)
1001         Dtilde(1,2,i)=b(8,i)
1002         Dtilde(2,1,i)=-b(8,i)
1003         Dtilde(2,2,i)=b(6,i)
1004         Dtilde(1,1,-i)=b(6,i)
1005         Dtilde(1,2,-i)=-b(8,i)
1006         Dtilde(2,1,-i)=b(8,i)
1007         Dtilde(2,2,-i)=b(6,i)
1008
1009 c        Dtilde(1,1,i)=0.0d0
1010 c        Dtilde(1,2,i)=0.0d0
1011 c        Dtilde(2,1,i)=0.0d0
1012 c        Dtilde(2,2,i)=0.0d0
1013         EEold(1,1,i)= b(10,i)+b(11,i)
1014         EEold(2,2,i)=-b(10,i)+b(11,i)
1015         EEold(2,1,i)= b(12,i)-b(13,i)
1016         EEold(1,2,i)= b(12,i)+b(13,i)
1017         EEold(1,1,-i)= b(10,i)+b(11,i)
1018         EEold(2,2,-i)=-b(10,i)+b(11,i)
1019         EEold(2,1,-i)=-b(12,i)+b(13,i)
1020         EEold(1,2,-i)=-b(12,i)-b(13,i)
1021         write(iout,*) "TU DOCHODZE"
1022         print *,"JESTEM"
1023 c        ee(1,1,i)=1.0d0
1024 c        ee(2,2,i)=1.0d0
1025 c        ee(2,1,i)=0.0d0
1026 c        ee(1,2,i)=0.0d0
1027 c        ee(2,1,i)=ee(1,2,i)
1028       enddo
1029 c      lprint=.true.
1030       if (lprint) then
1031       do i=1,nloctyp
1032         write (iout,*) 'Type',i
1033         write (iout,*) 'B1'
1034         write(iout,*) B1(1,i),B1(2,i)
1035         write (iout,*) 'B2'
1036         write(iout,*) B2(1,i),B2(2,i)
1037         write (iout,*) 'CC'
1038         do j=1,2
1039           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1040         enddo
1041         write(iout,*) 'DD'
1042         do j=1,2
1043           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1044         enddo
1045         write(iout,*) 'EE'
1046         do j=1,2
1047           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1048         enddo
1049       enddo
1050       endif
1051 c      lprint=.false.
1052
1053
1054 C Read electrostatic-interaction parameters
1055 C
1056       if (lprint) then
1057         write (iout,*)
1058         write (iout,'(/a)') 'Electrostatic interaction constants:'
1059         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
1060      &            'IT','JT','APP','BPP','AEL6','AEL3'
1061       endif
1062       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1063       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1064       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1065       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1066       close (ielep)
1067       do i=1,2
1068         do j=1,2
1069         rri=rpp(i,j)**6
1070         app (i,j)=epp(i,j)*rri*rri 
1071         bpp (i,j)=-2.0D0*epp(i,j)*rri
1072         ael6(i,j)=elpp6(i,j)*4.2D0**6
1073         ael3(i,j)=elpp3(i,j)*4.2D0**3
1074 c        lprint=.true.
1075         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1076      &                    ael6(i,j),ael3(i,j)
1077 c        lprint=.false.
1078         enddo
1079       enddo
1080 C
1081 C Read side-chain interaction parameters.
1082 C
1083       read (isidep,*,end=117,err=117) ipot,expon
1084       if (ipot.lt.1 .or. ipot.gt.5) then
1085         write (iout,'(2a)') 'Error while reading SC interaction',
1086      &               'potential file - unknown potential type.'
1087 #ifdef MPI
1088         call MPI_Finalize(Ierror)
1089 #endif
1090         stop
1091       endif
1092       expon2=expon/2
1093       if(me.eq.king)
1094      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1095      & ', exponents are ',expon,2*expon 
1096       goto (10,20,30,30,40) ipot
1097 C----------------------- LJ potential ---------------------------------
1098    10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1099      &   (sigma0(i),i=1,ntyp)
1100       if (lprint) then
1101         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1102         write (iout,'(a/)') 'The epsilon array:'
1103         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1104         write (iout,'(/a)') 'One-body parameters:'
1105         write (iout,'(a,4x,a)') 'residue','sigma'
1106         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1107       endif
1108       goto 50
1109 C----------------------- LJK potential --------------------------------
1110    20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1111      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1112       if (lprint) then
1113         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1114         write (iout,'(a/)') 'The epsilon array:'
1115         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1116         write (iout,'(/a)') 'One-body parameters:'
1117         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1118         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1119      &        i=1,ntyp)
1120       endif
1121       goto 50
1122 C---------------------- GB or BP potential -----------------------------
1123    30 do i=1,ntyp
1124        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1125       enddo
1126       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1127       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1128       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1129       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1130 C now we start reading lipid
1131       do i=1,ntyp
1132        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1133        print *,"WARNING!!"
1134        do j=1,ntyp
1135        epslip(i,j)=epslip(i,j)+0.05d0
1136        enddo
1137       enddo
1138 C For the GB potential convert sigma'**2 into chi'
1139       if (ipot.eq.4) then
1140         do i=1,ntyp
1141           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1142         enddo
1143       endif
1144       if (lprint) then
1145         write (iout,'(/a/)') 'Parameters of the BP potential:'
1146         write (iout,'(a/)') 'The epsilon array:'
1147         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1148         write (iout,'(/a)') 'One-body parameters:'
1149         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1150      &       '    chip  ','    alph  '
1151         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1152      &                     chip(i),alp(i),i=1,ntyp)
1153       endif
1154       goto 50
1155 C--------------------- GBV potential -----------------------------------
1156    40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1157      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1158      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1159       if (lprint) then
1160         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1161         write (iout,'(a/)') 'The epsilon array:'
1162         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1163         write (iout,'(/a)') 'One-body parameters:'
1164         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1165      &      's||/s_|_^2','    chip  ','    alph  '
1166         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1167      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1168       endif
1169    50 continue
1170       close (isidep)
1171 C-----------------------------------------------------------------------
1172 C Calculate the "working" parameters of SC interactions.
1173       do i=2,ntyp
1174         do j=1,i-1
1175           eps(i,j)=eps(j,i)
1176           epslip(i,j)=epslip(j,i)
1177         enddo
1178       enddo
1179       do i=1,ntyp
1180         do j=i,ntyp
1181           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1182           sigma(j,i)=sigma(i,j)
1183           rs0(i,j)=dwa16*sigma(i,j)
1184           rs0(j,i)=rs0(i,j)
1185         enddo
1186       enddo
1187       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
1188      & 'Working parameters of the SC interactions:',
1189      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1190      & '  chi1   ','   chi2   ' 
1191       do i=1,ntyp
1192         do j=i,ntyp
1193           epsij=eps(i,j)
1194           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1195             rrij=sigma(i,j)
1196           else
1197             rrij=rr0(i)+rr0(j)
1198           endif
1199           r0(i,j)=rrij
1200           r0(j,i)=rrij
1201           rrij=rrij**expon
1202           epsij=eps(i,j)
1203           sigeps=dsign(1.0D0,epsij)
1204           epsij=dabs(epsij)
1205           aa_aq(i,j)=epsij*rrij*rrij
1206           bb_aq(i,j)=-sigeps*epsij*rrij
1207           aa_aq(j,i)=aa_aq(i,j)
1208           bb_aq(j,i)=bb_aq(i,j)
1209           epsijlip=epslip(i,j)
1210           sigeps=dsign(1.0D0,epsijlip)
1211           epsijlip=dabs(epsijlip)
1212           aa_lip(i,j)=epsijlip*rrij*rrij
1213           bb_lip(i,j)=-sigeps*epsijlip*rrij
1214           aa_lip(j,i)=aa_lip(i,j)
1215           bb_lip(j,i)=bb_lip(i,j)
1216           if (ipot.gt.2) then
1217             sigt1sq=sigma0(i)**2
1218             sigt2sq=sigma0(j)**2
1219             sigii1=sigii(i)
1220             sigii2=sigii(j)
1221             ratsig1=sigt2sq/sigt1sq
1222             ratsig2=1.0D0/ratsig1
1223             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1224             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1225             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1226           else
1227             rsum_max=sigma(i,j)
1228           endif
1229 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1230             sigmaii(i,j)=rsum_max
1231             sigmaii(j,i)=rsum_max 
1232 c         else
1233 c           sigmaii(i,j)=r0(i,j)
1234 c           sigmaii(j,i)=r0(i,j)
1235 c         endif
1236 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1237           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1238             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1239             augm(i,j)=epsij*r_augm**(2*expon)
1240 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1241             augm(j,i)=augm(i,j)
1242           else
1243             augm(i,j)=0.0D0
1244             augm(j,i)=0.0D0
1245           endif
1246           if (lprint) then
1247             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
1248      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1249      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1250           endif
1251         enddo
1252       enddo
1253 #ifdef OLDSCP
1254 C
1255 C Define the SC-p interaction constants (hard-coded; old style)
1256 C
1257       do i=1,ntyp
1258 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1259 C helix formation)
1260 c       aad(i,1)=0.3D0*4.0D0**12
1261 C Following line for constants currently implemented
1262 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1263         aad(i,1)=1.5D0*4.0D0**12
1264 c       aad(i,1)=0.17D0*5.6D0**12
1265         aad(i,2)=aad(i,1)
1266 C "Soft" SC-p repulsion
1267         bad(i,1)=0.0D0
1268 C Following line for constants currently implemented
1269 c       aad(i,1)=0.3D0*4.0D0**6
1270 C "Hard" SC-p repulsion
1271         bad(i,1)=3.0D0*4.0D0**6
1272 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1273         bad(i,2)=bad(i,1)
1274 c       aad(i,1)=0.0D0
1275 c       aad(i,2)=0.0D0
1276 c       bad(i,1)=1228.8D0
1277 c       bad(i,2)=1228.8D0
1278       enddo
1279 #else
1280 C
1281 C 8/9/01 Read the SC-p interaction constants from file
1282 C
1283       do i=1,ntyp
1284         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1285       enddo
1286       do i=1,ntyp
1287         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1288         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1289         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1290         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1291       enddo
1292 c      lprint=.true.
1293       if (lprint) then
1294         write (iout,*) "Parameters of SC-p interactions:"
1295         do i=1,ntyp
1296           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1297      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1298         enddo
1299       endif
1300 c      lprint=.false.
1301 #endif
1302 C
1303 C Define the constants of the disulfide bridge
1304 C
1305 C      ebr=-12.00D0
1306 c
1307 c Old arbitrary potential - commented out.
1308 c
1309 c      dbr= 4.20D0
1310 c      fbr= 3.30D0
1311 c
1312 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1313 c energy surface of diethyl disulfide.
1314 c A. Liwo and U. Kozlowska, 11/24/03
1315 c
1316 C      D0CM = 3.78d0
1317 C      AKCM = 15.1d0
1318 C      AKTH = 11.0d0
1319 C      AKCT = 12.0d0
1320 C      V1SS =-1.08d0
1321 C      V2SS = 7.61d0
1322 C      V3SS = 13.7d0
1323 c      akcm=0.0d0
1324 c      akth=0.0d0
1325 c      akct=0.0d0
1326 c      v1ss=0.0d0
1327 c      v2ss=0.0d0
1328 c      v3ss=0.0d0
1329       
1330 C      if(me.eq.king) then
1331 C      write (iout,'(/a)') "Disulfide bridge parameters:"
1332 C      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1333 C      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1334 C      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1335 C      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1336 C     &  ' v3ss:',v3ss
1337 C      endif
1338 C set the variables used for shielding effect
1339 C      write (iout,*) "SHIELD MODE",shield_mode
1340 C      if (shield_mode.gt.0) then
1341 C VSolvSphere the volume of solving sphere
1342 C      print *,pi,"pi"
1343 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
1344 C there will be no distinction between proline peptide group and normal peptide
1345 C group in case of shielding parameters
1346 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1347 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1348 C      write (iout,*) VSolvSphere,VSolvSphere_div
1349 C long axis of side chain 
1350 C      do i=1,ntyp
1351 C      long_r_sidechain(i)=vbldsc0(1,i)
1352 C      short_r_sidechain(i)=sigma0(i)
1353 C      enddo
1354 C lets set the buffor value
1355 C      buff_shield=1.0d0
1356 C      endif
1357       return
1358   111 write (iout,*) "Error reading bending energy parameters."
1359       goto 999
1360   112 write (iout,*) "Error reading rotamer energy parameters."
1361       goto 999
1362   113 write (iout,*) "Error reading torsional energy parameters."
1363       goto 999
1364   114 write (iout,*) "Error reading double torsional energy parameters."
1365       goto 999
1366   115 write (iout,*) 
1367      &  "Error reading cumulant (multibody energy) parameters."
1368       goto 999
1369   116 write (iout,*) "Error reading electrostatic energy parameters."
1370       goto 999
1371  1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1372       goto 999
1373   117 write (iout,*) "Error reading side chain interaction parameters."
1374       goto 999
1375   118 write (iout,*) "Error reading SCp interaction parameters."
1376       goto 999
1377   119 write (iout,*) "Error reading SCCOR parameters"
1378   999 continue
1379 #ifdef MPI
1380       call MPI_Finalize(Ierror)
1381 #endif
1382       stop
1383       return
1384       end
1385
1386
1387       subroutine getenv_loc(var, val)
1388       character(*) var, val
1389
1390 #ifdef WINIFL
1391       character(2000) line
1392       external ilen
1393
1394       open (196,file='env',status='old',readonly,shared)
1395       iread=0
1396 c      write(*,*)'looking for ',var
1397 10    read(196,*,err=11,end=11)line
1398       iread=index(line,var)
1399 c      write(*,*)iread,' ',var,' ',line
1400       if (iread.eq.0) go to 10 
1401 c      write(*,*)'---> ',line
1402 11    continue
1403       if(iread.eq.0) then
1404 c       write(*,*)'CHUJ'
1405        val=''
1406       else
1407        iread=iread+ilen(var)+1
1408        read (line(iread:),*,err=12,end=12) val
1409 c       write(*,*)'OK: ',var,' = ',val
1410       endif
1411       close(196)
1412       return
1413 12    val=''
1414       close(196)
1415 #elif (defined CRAY)
1416       integer lennam,lenval,ierror
1417 c
1418 c        getenv using a POSIX call, useful on the T3D
1419 c        Sept 1996, comment out error check on advice of H. Pritchard
1420 c
1421       lennam = len(var)
1422       if(lennam.le.0) stop '--error calling getenv--'
1423       call pxfgetenv(var,lennam,val,lenval,ierror)
1424 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1425 #else
1426       call getenv(var,val)
1427 #endif
1428 C set the variables used for shielding effect
1429 C      if (shield_mode.gt.0) then
1430 C VSolvSphere the volume of solving sphere
1431 C      print *,pi,"pi"
1432 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
1433 C there will be no distinction between proline peptide group and normal peptide
1434 C group in case of shielding parameters
1435 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1436 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1437 C long axis of side chain 
1438 C      do i=1,ntyp
1439 C      long_r_sidechain(i)=vbldsc0(1,i)
1440 C      short_r_sidechain(i)=sigma0(i)
1441 C      enddo
1442 C lets set the buffor value
1443 C      buff_shield=1.0d0
1444 C      endif
1445       return
1446       end