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