arcos = 0.5D0*(PI-DSIGN(1.0D0,X)*PI)
[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,"SCCORPAR",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 Read of Side-chain backbone correlation parameters
586 C Modified 11 May 2012 by Adasko
587 CCC
588 C
589       read (isccor,*) nsccortyp
590       read (isccor,*) (isccortyp(i),i=1,ntyp)
591       do i=-ntyp,-1
592         isccortyp(i)=-isccortyp(-i)
593       enddo
594       iscprol=isccortyp(20)
595 c      write (iout,*) 'ntortyp',ntortyp
596       maxinter=3
597 cc maxinter is maximum interaction sites
598       do l=1,maxinter
599       do i=1,nsccortyp
600         do j=1,nsccortyp
601           read (isccor,*)
602      &nterm_sccor(i,j),nlor_sccor(i,j)
603           write (iout,*) nterm_sccor(i,j)
604           v0ijsccor=0.0d0
605           v0ijsccor1=0.0d0
606           v0ijsccor2=0.0d0
607           v0ijsccor3=0.0d0
608           si=-1.0d0
609           nterm_sccor(-i,j)=nterm_sccor(i,j)
610           nterm_sccor(-i,-j)=nterm_sccor(i,j)
611           nterm_sccor(i,-j)=nterm_sccor(i,j)
612           write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
613      &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
614           do k=1,nterm_sccor(i,j)
615             read (isccor,*) kk,v1sccor(k,l,i,j)
616      &    ,v2sccor(k,l,i,j)
617             if (j.eq.iscprol) then
618              if (i.eq.isccortyp(10)) then
619              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
620              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
621              else
622              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
623      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
624              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
625      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
626              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
627              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
628              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
629              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
630              endif
631             else
632              if (i.eq.isccortyp(10)) then
633              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
634              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
635              else
636                if (j.eq.isccortyp(10)) then
637              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
638              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
639                else
640              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
641              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
642              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
643              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
644              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
645              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
646                 endif
647                endif
648             endif
649             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
650             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
651             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
652             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
653             si=-si
654            enddo
655           do k=1,nlor_sccor(i,j)
656             read (isccor,*) kk,vlor1sccor(k,i,j),
657      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
658             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
659      &(1+vlor3sccor(k,i,j)**2)
660           enddo
661           v0sccor(l,i,j)=v0ijsccor
662           v0sccor(l,-i,j)=v0ijsccor1
663           v0sccor(l,i,-j)=v0ijsccor2
664           v0sccor(l,-i,-j)=v0ijsccor3
665           enddo
666         enddo
667       enddo
668       close (isccor)
669       if (lprint) then
670         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
671         do i=1,nsccortyp
672           do j=1,nsccortyp
673             write (iout,*) 'ityp',i,' jtyp',j
674             write (iout,*) 'Fourier constants'
675             do k=1,nterm_sccor(i,j)
676               write (iout,'(2(1pe15.5))')
677      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
678             enddo
679             write (iout,*) 'Lorenz constants'
680             do k=1,nlor_sccor(i,j)
681               write (iout,'(3(1pe15.5))')
682      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
683             enddo
684           enddo
685         enddo
686       endif
687 C
688 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
689 C         interaction energy of the Gly, Ala, and Pro prototypes.
690 C
691       read (ifourier,*) nloctyp
692       do i=1,nloctyp
693         read (ifourier,*)
694         read (ifourier,*) (b(ii,i),ii=1,13)
695         if (lprint) then
696         write (iout,*) 'Type',i
697         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
698         endif
699         B1(1,i)  = b(3,i)
700         B1(2,i)  = b(5,i)
701         B1tilde(1,i) = b(3,i)
702         B1tilde(2,i) =-b(5,i) 
703         B2(1,i)  = b(2,i)
704         B2(2,i)  = b(4,i)
705         CC(1,1,i)= b(7,i)
706         CC(2,2,i)=-b(7,i)
707         CC(2,1,i)= b(9,i)
708         CC(1,2,i)= b(9,i)
709         Ctilde(1,1,i)=b(7,i)
710         Ctilde(1,2,i)=b(9,i)
711         Ctilde(2,1,i)=-b(9,i)
712         Ctilde(2,2,i)=b(7,i)
713         DD(1,1,i)= b(6,i)
714         DD(2,2,i)=-b(6,i)
715         DD(2,1,i)= b(8,i)
716         DD(1,2,i)= b(8,i)
717         Dtilde(1,1,i)=b(6,i)
718         Dtilde(1,2,i)=b(8,i)
719         Dtilde(2,1,i)=-b(8,i)
720         Dtilde(2,2,i)=b(6,i)
721         EE(1,1,i)= b(10,i)+b(11,i)
722         EE(2,2,i)=-b(10,i)+b(11,i)
723         EE(2,1,i)= b(12,i)-b(13,i)
724         EE(1,2,i)= b(12,i)+b(13,i)
725       enddo
726       if (lprint) then
727       do i=1,nloctyp
728         write (iout,*) 'Type',i
729         write (iout,*) 'B1'
730 c        write (iout,'(f10.5)') B1(:,i)
731         write(iout,*) B1(1,i),B1(2,i)
732         write (iout,*) 'B2'
733 c        write (iout,'(f10.5)') B2(:,i)
734         write(iout,*) B2(1,i),B2(2,i)
735         write (iout,*) 'CC'
736         do j=1,2
737           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
738         enddo
739         write(iout,*) 'DD'
740         do j=1,2
741           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
742         enddo
743         write(iout,*) 'EE'
744         do j=1,2
745           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
746         enddo
747       enddo
748       endif
749
750 C Read electrostatic-interaction parameters
751 C
752       if (lprint) then
753         write (iout,'(/a)') 'Electrostatic interaction constants:'
754         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
755      &            'IT','JT','APP','BPP','AEL6','AEL3'
756       endif
757       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
758       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
759       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
760       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
761       close (ielep)
762       do i=1,2
763         do j=1,2
764         rri=rpp(i,j)**6
765         app (i,j)=epp(i,j)*rri*rri 
766         bpp (i,j)=-2.0D0*epp(i,j)*rri
767         ael6(i,j)=elpp6(i,j)*4.2D0**6
768         ael3(i,j)=elpp3(i,j)*4.2D0**3
769         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
770      &                    ael6(i,j),ael3(i,j)
771         enddo
772       enddo
773 C
774 C Read side-chain interaction parameters.
775 C
776       read (isidep,*) ipot,expon
777       if (ipot.lt.1 .or. ipot.gt.5) then
778         write (iout,'(2a)') 'Error while reading SC interaction',
779      &               'potential file - unknown potential type.'
780         stop
781       endif
782       expon2=expon/2
783       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
784      & ', exponents are ',expon,2*expon 
785       goto (10,20,30,30,40) ipot
786 C----------------------- LJ potential ---------------------------------
787    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
788       if (lprint) then
789         write (iout,'(/a/)') 'Parameters of the LJ potential:'
790         write (iout,'(a/)') 'The epsilon array:'
791         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
792         write (iout,'(/a)') 'One-body parameters:'
793         write (iout,'(a,4x,a)') 'residue','sigma'
794         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
795       endif
796       goto 50
797 C----------------------- LJK potential --------------------------------
798    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
799      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
800       if (lprint) then
801         write (iout,'(/a/)') 'Parameters of the LJK potential:'
802         write (iout,'(a/)') 'The epsilon array:'
803         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
804         write (iout,'(/a)') 'One-body parameters:'
805         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
806         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
807      &        i=1,ntyp)
808       endif
809       goto 50
810 C---------------------- GB or BP potential -----------------------------
811    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
812      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
813      &  (alp(i),i=1,ntyp)
814 C For the GB potential convert sigma'**2 into chi'
815       if (ipot.eq.4) then
816         do i=1,ntyp
817           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
818         enddo
819       endif
820       if (lprint) then
821         write (iout,'(/a/)') 'Parameters of the BP potential:'
822         write (iout,'(a/)') 'The epsilon array:'
823         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
824         write (iout,'(/a)') 'One-body parameters:'
825         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
826      &       '    chip  ','    alph  '
827         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
828      &                     chip(i),alp(i),i=1,ntyp)
829       endif
830       goto 50
831 C--------------------- GBV potential -----------------------------------
832    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
833      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
834      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
835       if (lprint) then
836         write (iout,'(/a/)') 'Parameters of the GBV potential:'
837         write (iout,'(a/)') 'The epsilon array:'
838         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
839         write (iout,'(/a)') 'One-body parameters:'
840         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
841      &      's||/s_|_^2','    chip  ','    alph  '
842         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
843      &           sigii(i),chip(i),alp(i),i=1,ntyp)
844       endif
845    50 continue
846       close (isidep)
847 C-----------------------------------------------------------------------
848 C Calculate the "working" parameters of SC interactions.
849       do i=2,ntyp
850         do j=1,i-1
851           eps(i,j)=eps(j,i)
852         enddo
853       enddo
854       do i=1,ntyp
855         do j=i,ntyp
856           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
857           sigma(j,i)=sigma(i,j)
858           rs0(i,j)=dwa16*sigma(i,j)
859           rs0(j,i)=rs0(i,j)
860         enddo
861       enddo
862       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
863      & 'Working parameters of the SC interactions:',
864      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
865      & '  chi1   ','   chi2   ' 
866       do i=1,ntyp
867         do j=i,ntyp
868           epsij=eps(i,j)
869           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
870             rrij=sigma(i,j)
871           else
872             rrij=rr0(i)+rr0(j)
873           endif
874           r0(i,j)=rrij
875           r0(j,i)=rrij
876           rrij=rrij**expon
877           epsij=eps(i,j)
878           sigeps=dsign(1.0D0,epsij)
879           epsij=dabs(epsij)
880           aa(i,j)=epsij*rrij*rrij
881           bb(i,j)=-sigeps*epsij*rrij
882           aa(j,i)=aa(i,j)
883           bb(j,i)=bb(i,j)
884           if (ipot.gt.2) then
885             sigt1sq=sigma0(i)**2
886             sigt2sq=sigma0(j)**2
887             sigii1=sigii(i)
888             sigii2=sigii(j)
889             ratsig1=sigt2sq/sigt1sq
890             ratsig2=1.0D0/ratsig1
891             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
892             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
893             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
894           else
895             rsum_max=sigma(i,j)
896           endif
897 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
898             sigmaii(i,j)=rsum_max
899             sigmaii(j,i)=rsum_max 
900 c         else
901 c           sigmaii(i,j)=r0(i,j)
902 c           sigmaii(j,i)=r0(i,j)
903 c         endif
904 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
905           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
906             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
907             augm(i,j)=epsij*r_augm**(2*expon)
908 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
909             augm(j,i)=augm(i,j)
910           else
911             augm(i,j)=0.0D0
912             augm(j,i)=0.0D0
913           endif
914           if (lprint) then
915             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
916      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
917      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
918           endif
919         enddo
920       enddo
921 C
922 C Define the SC-p interaction constants
923 C
924 #ifdef OLDSCP
925       do i=1,20
926 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
927 C helix formation)
928 c       aad(i,1)=0.3D0*4.0D0**12
929 C Following line for constants currently implemented
930 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
931         aad(i,1)=1.5D0*4.0D0**12
932 c       aad(i,1)=0.17D0*5.6D0**12
933         aad(i,2)=aad(i,1)
934 C "Soft" SC-p repulsion
935         bad(i,1)=0.0D0
936 C Following line for constants currently implemented
937 c       aad(i,1)=0.3D0*4.0D0**6
938 C "Hard" SC-p repulsion
939         bad(i,1)=3.0D0*4.0D0**6
940 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
941         bad(i,2)=bad(i,1)
942 c       aad(i,1)=0.0D0
943 c       aad(i,2)=0.0D0
944 c       bad(i,1)=1228.8D0
945 c       bad(i,2)=1228.8D0
946       enddo
947 #else
948 C
949 C 8/9/01 Read the SC-p interaction constants from file
950 C
951       do i=1,ntyp
952         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
953       enddo
954       do i=1,ntyp
955         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
956         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
957         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
958         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
959       enddo
960
961       if (lprint) then
962         write (iout,*) "Parameters of SC-p interactions:"
963         do i=1,20
964           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
965      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
966         enddo
967       endif
968 #endif
969 C
970 C Define the constants of the disulfide bridge
971 C
972       ebr=-5.50D0
973 c
974 c Old arbitrary potential - commented out.
975 c
976 c      dbr= 4.20D0
977 c      fbr= 3.30D0
978 c
979 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
980 c energy surface of diethyl disulfide.
981 c A. Liwo and U. Kozlowska, 11/24/03
982 c
983       D0CM = 3.78d0
984       AKCM = 15.1d0
985       AKTH = 11.0d0
986       AKCT = 12.0d0
987       V1SS =-1.08d0
988       V2SS = 7.61d0
989       V3SS = 13.7d0
990
991       if (lprint) then
992       write (iout,'(/a)') "Disulfide bridge parameters:"
993       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
994       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
995       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
996       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
997      & ' v3ss:',v3ss
998       endif
999       return
1000       end