8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2,ijunk
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(ntyp+1,5)) !(ntyp+1)
817 allocate(isc(ntyp+1,5)) !(ntyp+1)
818 allocate(restok(ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
871 read(iion,*) msc(i,5),restok(i,5)
872 print *,msc(i,5),restok(i,5)
876 read (iion,*) ncatprotparm
877 allocate(catprm(ncatprotparm,4))
879 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
882 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
883 !----------------------------------------------------
884 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
885 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
886 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
887 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
888 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
889 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
899 allocate(liptranene(ntyp))
900 !C reading lipid parameters
901 write (iout,*) "iliptranpar",iliptranpar
903 read(iliptranpar,*) pepliptran
906 read(iliptranpar,*) liptranene(i)
907 print *,liptranene(i)
913 ! Read the parameters of the probability distribution/energy expression
914 ! of the virtual-bond valence angles theta
917 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
918 (bthet(j,i,1,1),j=1,2)
919 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
920 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
921 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
925 athet(1,i,1,-1)=athet(1,i,1,1)
926 athet(2,i,1,-1)=athet(2,i,1,1)
927 bthet(1,i,1,-1)=-bthet(1,i,1,1)
928 bthet(2,i,1,-1)=-bthet(2,i,1,1)
929 athet(1,i,-1,1)=-athet(1,i,1,1)
930 athet(2,i,-1,1)=-athet(2,i,1,1)
931 bthet(1,i,-1,1)=bthet(1,i,1,1)
932 bthet(2,i,-1,1)=bthet(2,i,1,1)
936 athet(1,i,-1,-1)=athet(1,-i,1,1)
937 athet(2,i,-1,-1)=-athet(2,-i,1,1)
938 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
939 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
940 athet(1,i,-1,1)=athet(1,-i,1,1)
941 athet(2,i,-1,1)=-athet(2,-i,1,1)
942 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
943 bthet(2,i,-1,1)=bthet(2,-i,1,1)
944 athet(1,i,1,-1)=-athet(1,-i,1,1)
945 athet(2,i,1,-1)=athet(2,-i,1,1)
946 bthet(1,i,1,-1)=bthet(1,-i,1,1)
947 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
952 polthet(j,i)=polthet(j,-i)
955 gthet(j,i)=gthet(j,-i)
963 'Parameters of the virtual-bond valence angles:'
964 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
965 ' ATHETA0 ',' A1 ',' A2 ',&
968 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
969 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
971 write (iout,'(/a/9x,5a/79(1h-))') &
972 'Parameters of the expression for sigma(theta_c):',&
973 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
974 ' ALPH3 ',' SIGMA0C '
976 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
977 (polthet(j,i),j=0,3),sigc0(i)
979 write (iout,'(/a/9x,5a/79(1h-))') &
980 'Parameters of the second gaussian:',&
981 ' THETA0 ',' SIGMA0 ',' G1 ',&
984 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
985 sig0(i),(gthet(j,i),j=1,3)
989 'Parameters of the virtual-bond valence angles:'
990 write (iout,'(/a/9x,5a/79(1h-))') &
991 'Coefficients of expansion',&
992 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
993 ' b1*10^1 ',' b2*10^1 '
995 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
996 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
997 (10*bthet(j,i,1,1),j=1,2)
999 write (iout,'(/a/9x,5a/79(1h-))') &
1000 'Parameters of the expression for sigma(theta_c):',&
1001 ' alpha0 ',' alph1 ',' alph2 ',&
1002 ' alhp3 ',' sigma0c '
1004 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1005 (polthet(j,i),j=0,3),sigc0(i)
1007 write (iout,'(/a/9x,5a/79(1h-))') &
1008 'Parameters of the second gaussian:',&
1009 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1012 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1013 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1019 ! Read the parameters of Utheta determined from ab initio surfaces
1020 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1022 IF (tor_mode.eq.0) THEN
1023 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1024 ntheterm3,nsingle,ndouble
1025 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1027 !----------------------------------------------------
1028 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1029 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1030 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1031 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1032 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1033 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1034 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1035 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1036 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1037 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1038 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1039 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1040 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1041 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1042 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1043 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1044 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1045 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1046 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1047 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1048 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1049 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1050 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1051 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1053 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1055 ithetyp(i)=-ithetyp(-i)
1058 aa0thet(:,:,:,:)=0.0d0
1059 aathet(:,:,:,:,:)=0.0d0
1060 bbthet(:,:,:,:,:,:)=0.0d0
1061 ccthet(:,:,:,:,:,:)=0.0d0
1062 ddthet(:,:,:,:,:,:)=0.0d0
1063 eethet(:,:,:,:,:,:)=0.0d0
1064 ffthet(:,:,:,:,:,:,:)=0.0d0
1065 ggthet(:,:,:,:,:,:,:)=0.0d0
1067 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1069 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1070 ! VAR:1=non-glicyne non-proline 2=proline
1071 ! VAR:negative values for D-aminoacid
1073 do j=-nthetyp,nthetyp
1074 do k=-nthetyp,nthetyp
1075 read (ithep,'(6a)',end=111,err=111) res1
1076 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1077 ! VAR: aa0thet is variable describing the average value of Foureir
1078 ! VAR: expansion series
1079 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1080 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1081 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1082 read (ithep,*,end=111,err=111) &
1083 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1084 read (ithep,*,end=111,err=111) &
1085 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1086 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1087 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1088 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1090 read (ithep,*,end=111,err=111) &
1091 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1092 ffthet(lll,llll,ll,i,j,k,iblock),&
1093 ggthet(llll,lll,ll,i,j,k,iblock),&
1094 ggthet(lll,llll,ll,i,j,k,iblock),&
1095 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1100 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1101 ! coefficients of theta-and-gamma-dependent terms are zero.
1102 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1103 ! RECOMENTDED AFTER VERSION 3.3)
1107 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1108 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1110 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1111 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1114 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1116 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1119 ! AND COMMENT THE LOOPS BELOW
1123 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1124 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1126 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1127 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1130 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1132 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1137 ! Substitution for D aminoacids from symmetry.
1140 do j=-nthetyp,nthetyp
1141 do k=-nthetyp,nthetyp
1142 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1144 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1148 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1149 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1150 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1151 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1157 ffthet(llll,lll,ll,i,j,k,iblock)= &
1158 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1159 ffthet(lll,llll,ll,i,j,k,iblock)= &
1160 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1161 ggthet(llll,lll,ll,i,j,k,iblock)= &
1162 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1163 ggthet(lll,llll,ll,i,j,k,iblock)= &
1164 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1173 ! Control printout of the coefficients of virtual-bond-angle potentials
1176 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1181 write (iout,'(//4a)') &
1182 'Type ',onelett(i),onelett(j),onelett(k)
1183 write (iout,'(//a,10x,a)') " l","a[l]"
1184 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1185 write (iout,'(i2,1pe15.5)') &
1186 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1188 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1189 "b",l,"c",l,"d",l,"e",l
1191 write (iout,'(i2,4(1pe15.5))') m,&
1192 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1193 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1197 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1198 "f+",l,"f-",l,"g+",l,"g-",l
1201 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1202 ffthet(n,m,l,i,j,k,iblock),&
1203 ffthet(m,n,l,i,j,k,iblock),&
1204 ggthet(n,m,l,i,j,k,iblock),&
1205 ggthet(m,n,l,i,j,k,iblock)
1216 !C here will be the apropriate recalibrating for D-aminoacid
1217 read (ithep,*,end=121,err=121) nthetyp
1218 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1219 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1220 do i=-nthetyp+1,nthetyp-1
1221 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1222 do j=0,nbend_kcc_Tb(i)
1223 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1227 write (iout,'(a)') &
1228 "Parameters of the valence-only potentials"
1229 do i=-nthetyp+1,nthetyp-1
1230 write (iout,'(2a)') "Type ",toronelet(i)
1231 do k=0,nbend_kcc_Tb(i)
1232 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1238 write (2,*) "Start reading THETA_PDB",ithep_pdb
1240 ! write (2,*) 'i=',i
1241 read (ithep_pdb,*,err=111,end=111) &
1242 a0thet(i),(athet(j,i,1,1),j=1,2),&
1243 (bthet(j,i,1,1),j=1,2)
1244 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1245 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1246 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1247 sigc0(i)=sigc0(i)**2
1250 athet(1,i,1,-1)=athet(1,i,1,1)
1251 athet(2,i,1,-1)=athet(2,i,1,1)
1252 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1253 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1254 athet(1,i,-1,1)=-athet(1,i,1,1)
1255 athet(2,i,-1,1)=-athet(2,i,1,1)
1256 bthet(1,i,-1,1)=bthet(1,i,1,1)
1257 bthet(2,i,-1,1)=bthet(2,i,1,1)
1260 a0thet(i)=a0thet(-i)
1261 athet(1,i,-1,-1)=athet(1,-i,1,1)
1262 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1263 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1264 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1265 athet(1,i,-1,1)=athet(1,-i,1,1)
1266 athet(2,i,-1,1)=-athet(2,-i,1,1)
1267 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1268 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1269 athet(1,i,1,-1)=-athet(1,-i,1,1)
1270 athet(2,i,1,-1)=athet(2,-i,1,1)
1271 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1272 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1273 theta0(i)=theta0(-i)
1277 polthet(j,i)=polthet(j,-i)
1280 gthet(j,i)=gthet(j,-i)
1283 write (2,*) "End reading THETA_PDB"
1287 !--------------- Reading theta parameters for nucleic acid-------
1288 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1289 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1290 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1291 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1292 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1293 nthetyp_nucl+1,nthetyp_nucl+1))
1294 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1295 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1296 nthetyp_nucl+1,nthetyp_nucl+1))
1297 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1298 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1299 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1300 nthetyp_nucl+1,nthetyp_nucl+1))
1301 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1302 nthetyp_nucl+1,nthetyp_nucl+1))
1303 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1304 nthetyp_nucl+1,nthetyp_nucl+1))
1305 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1306 nthetyp_nucl+1,nthetyp_nucl+1))
1307 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1308 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1309 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1310 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1311 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1312 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1314 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1315 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1317 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1319 aa0thet_nucl(:,:,:)=0.0d0
1320 aathet_nucl(:,:,:,:)=0.0d0
1321 bbthet_nucl(:,:,:,:,:)=0.0d0
1322 ccthet_nucl(:,:,:,:,:)=0.0d0
1323 ddthet_nucl(:,:,:,:,:)=0.0d0
1324 eethet_nucl(:,:,:,:,:)=0.0d0
1325 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1326 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1331 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1332 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1333 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1334 read (ithep_nucl,*,end=111,err=111) &
1335 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1336 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1337 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1338 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1339 read (ithep_nucl,*,end=111,err=111) &
1340 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1341 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1342 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1347 !-------------------------------------------
1348 allocate(nlob(ntyp1)) !(ntyp1)
1349 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1350 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1351 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1358 gaussc(:,:,:,:)=0.0D0
1362 ! Read the parameters of the probability distribution/energy expression
1363 ! of the side chains.
1366 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1370 dsc_inv(i)=1.0D0/dsc(i)
1381 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1382 ((blower(k,l,1),l=1,k),k=1,3)
1383 censc(1,1,-i)=censc(1,1,i)
1384 censc(2,1,-i)=censc(2,1,i)
1385 censc(3,1,-i)=-censc(3,1,i)
1387 read (irotam,*,end=112,err=112) bsc(j,i)
1388 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1389 ((blower(k,l,j),l=1,k),k=1,3)
1390 censc(1,j,-i)=censc(1,j,i)
1391 censc(2,j,-i)=censc(2,j,i)
1392 censc(3,j,-i)=-censc(3,j,i)
1393 ! BSC is amplitude of Gaussian
1400 akl=akl+blower(k,m,j)*blower(l,m,j)
1404 if (((k.eq.3).and.(l.ne.3)) &
1405 .or.((l.eq.3).and.(k.ne.3))) then
1406 gaussc(k,l,j,-i)=-akl
1407 gaussc(l,k,j,-i)=-akl
1409 gaussc(k,l,j,-i)=akl
1410 gaussc(l,k,j,-i)=akl
1419 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1422 if (nlobi.gt.0) then
1424 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1425 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1426 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1427 'log h',(bsc(j,i),j=1,nlobi)
1428 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1429 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1431 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1432 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1435 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1436 write (iout,'(a,f10.4,4(16x,f10.4))') &
1437 'Center ',(bsc(j,i),j=1,nlobi)
1438 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1447 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1448 ! added by Urszula Kozlowska 07/11/2007
1450 !el Maximum number of SC local term fitting function coefficiants
1451 !el integer,parameter :: maxsccoef=65
1453 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1456 read (irotam,*,end=112,err=112)
1458 read (irotam,*,end=112,err=112)
1461 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1465 !---------reading nucleic acid parameters for rotamers-------------------
1466 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1467 do i=1,ntyp_molec(2)
1468 read (irotam_nucl,*,end=112,err=112)
1470 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1476 write (iout,*) "Base rotamer parameters"
1477 do i=1,ntyp_molec(2)
1478 write (iout,'(a)') restyp(i,2)
1479 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1484 ! Read the parameters of the probability distribution/energy expression
1485 ! of the side chains.
1487 write (2,*) "Start reading ROTAM_PDB"
1489 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1493 dsc_inv(i)=1.0D0/dsc(i)
1504 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1505 ((blower(k,l,1),l=1,k),k=1,3)
1507 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1508 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1509 ((blower(k,l,j),l=1,k),k=1,3)
1516 akl=akl+blower(k,m,j)*blower(l,m,j)
1526 write (2,*) "End reading ROTAM_PDB"
1532 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1533 !C interaction energy of the Gly, Ala, and Pro prototypes.
1535 read (ifourier,*) nloctyp
1536 SPLIT_FOURIERTOR = nloctyp.lt.0
1537 nloctyp = iabs(nloctyp)
1538 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1539 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1540 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1541 !C allocate(ctilde(2,2,nres))
1542 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1543 !C allocate(gtb1(2,nres))
1544 !C allocate(gtb2(2,nres))
1545 !C allocate(cc(2,2,nres))
1546 !C allocate(dd(2,2,nres))
1547 !C allocate(ee(2,2,nres))
1548 !C allocate(gtcc(2,2,nres))
1549 !C allocate(gtdd(2,2,nres))
1550 !C allocate(gtee(2,2,nres))
1553 allocate(itype2loc(-ntyp1:ntyp1))
1554 allocate(iloctyp(-nloctyp:nloctyp))
1555 allocate(bnew1(3,2,-nloctyp:nloctyp))
1556 allocate(bnew2(3,2,-nloctyp:nloctyp))
1557 allocate(ccnew(3,2,-nloctyp:nloctyp))
1558 allocate(ddnew(3,2,-nloctyp:nloctyp))
1559 allocate(e0new(3,-nloctyp:nloctyp))
1560 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1561 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1562 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1563 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1564 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1565 allocate(e0newtor(3,-nloctyp:nloctyp))
1566 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1568 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1569 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1570 itype2loc(ntyp1)=nloctyp
1571 iloctyp(nloctyp)=ntyp1
1573 itype2loc(-i)=-itype2loc(i)
1576 allocate(iloctyp(-nloctyp:nloctyp))
1583 iloctyp(-i)=-iloctyp(i)
1585 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1586 !c write (iout,*) "nloctyp",nloctyp,
1587 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1588 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1589 !c write (iout,*) "nloctyp",nloctyp,
1590 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1593 !c write (iout,*) "NEWCORR",i
1594 read (ifourier,*,end=115,err=115)
1597 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1600 !c write (iout,*) "NEWCORR BNEW1"
1601 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1604 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1607 !c write (iout,*) "NEWCORR BNEW2"
1608 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1610 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1611 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1613 !c write (iout,*) "NEWCORR CCNEW"
1614 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1616 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1617 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1619 !c write (iout,*) "NEWCORR DDNEW"
1620 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1624 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1628 !c write (iout,*) "NEWCORR EENEW1"
1629 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1631 read (ifourier,*,end=115,err=115) e0new(ii,i)
1633 !c write (iout,*) (e0new(ii,i),ii=1,3)
1635 !c write (iout,*) "NEWCORR EENEW"
1638 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1639 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1640 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1641 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1646 bnew1(ii,1,-i)= bnew1(ii,1,i)
1647 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1648 bnew2(ii,1,-i)= bnew2(ii,1,i)
1649 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1652 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1653 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1654 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1655 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1656 ccnew(ii,1,-i)=ccnew(ii,1,i)
1657 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1658 ddnew(ii,1,-i)=ddnew(ii,1,i)
1659 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1661 e0new(1,-i)= e0new(1,i)
1662 e0new(2,-i)=-e0new(2,i)
1663 e0new(3,-i)=-e0new(3,i)
1665 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1666 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1667 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1668 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1672 write (iout,'(a)') "Coefficients of the multibody terms"
1673 do i=-nloctyp+1,nloctyp-1
1674 write (iout,*) "Type: ",onelet(iloctyp(i))
1675 write (iout,*) "Coefficients of the expansion of B1"
1677 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1679 write (iout,*) "Coefficients of the expansion of B2"
1681 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1683 write (iout,*) "Coefficients of the expansion of C"
1684 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1685 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1686 write (iout,*) "Coefficients of the expansion of D"
1687 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1688 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1689 write (iout,*) "Coefficients of the expansion of E"
1690 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1693 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1698 IF (SPLIT_FOURIERTOR) THEN
1700 !c write (iout,*) "NEWCORR TOR",i
1701 read (ifourier,*,end=115,err=115)
1704 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1707 !c write (iout,*) "NEWCORR BNEW1 TOR"
1708 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1711 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1714 !c write (iout,*) "NEWCORR BNEW2 TOR"
1715 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1717 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1718 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1720 !c write (iout,*) "NEWCORR CCNEW TOR"
1721 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1723 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1724 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1726 !c write (iout,*) "NEWCORR DDNEW TOR"
1727 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1731 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1735 !c write (iout,*) "NEWCORR EENEW1 TOR"
1736 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1738 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1740 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1742 !c write (iout,*) "NEWCORR EENEW TOR"
1745 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1746 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1747 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1748 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1753 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1754 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1755 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1756 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1759 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1760 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1761 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1762 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1764 e0newtor(1,-i)= e0newtor(1,i)
1765 e0newtor(2,-i)=-e0newtor(2,i)
1766 e0newtor(3,-i)=-e0newtor(3,i)
1768 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1769 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1770 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1771 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1775 write (iout,'(a)') &
1776 "Single-body coefficients of the torsional potentials"
1777 do i=-nloctyp+1,nloctyp-1
1778 write (iout,*) "Type: ",onelet(iloctyp(i))
1779 write (iout,*) "Coefficients of the expansion of B1tor"
1781 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1782 j,(bnew1tor(k,j,i),k=1,3)
1784 write (iout,*) "Coefficients of the expansion of B2tor"
1786 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1787 j,(bnew2tor(k,j,i),k=1,3)
1789 write (iout,*) "Coefficients of the expansion of Ctor"
1790 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1791 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1792 write (iout,*) "Coefficients of the expansion of Dtor"
1793 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1794 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1795 write (iout,*) "Coefficients of the expansion of Etor"
1796 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1799 write (iout,'(1hE,2i1,2f10.5)') &
1800 j,k,(eenewtor(l,j,k,i),l=1,2)
1806 do i=-nloctyp+1,nloctyp-1
1809 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1814 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1818 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1819 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1820 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1821 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1824 ENDIF !SPLIT_FOURIER_TOR
1826 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1827 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1828 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1829 allocate(b(13,-nloctyp-1:nloctyp+1))
1831 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1833 read (ifourier,*,end=115,err=115)
1834 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1836 write (iout,*) 'Type ',onelet(iloctyp(i))
1837 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1851 !c B1tilde(1,i) = b(3)
1852 !c! B1tilde(2,i) =-b(5)
1853 !c! B1tilde(1,-i) =-b(3)
1854 !c! B1tilde(2,-i) =b(5)
1855 !c! b1tilde(1,i)=0.0d0
1856 !c b1tilde(2,i)=0.0d0
1861 !cc B1tilde(1,i) = b(3,i)
1862 !cc B1tilde(2,i) =-b(5,i)
1863 !c B1tilde(1,-i) =-b(3,i)
1864 !c B1tilde(2,-i) =b(5,i)
1865 !cc b1tilde(1,i)=0.0d0
1866 !cc b1tilde(2,i)=0.0d0
1867 !cc B2(1,i) = b(2,i)
1868 !cc B2(2,i) = b(4,i)
1870 !c B2(2,-i) =-b(4,i)
1874 CCold(1,1,i)= b(7,i)
1875 CCold(2,2,i)=-b(7,i)
1876 CCold(2,1,i)= b(9,i)
1877 CCold(1,2,i)= b(9,i)
1878 CCold(1,1,-i)= b(7,i)
1879 CCold(2,2,-i)=-b(7,i)
1880 CCold(2,1,-i)=-b(9,i)
1881 CCold(1,2,-i)=-b(9,i)
1886 !c Ctilde(1,1,i)= CCold(1,1,i)
1887 !c Ctilde(1,2,i)= CCold(1,2,i)
1888 !c Ctilde(2,1,i)=-CCold(2,1,i)
1889 !c Ctilde(2,2,i)=-CCold(2,2,i)
1894 !c Ctilde(1,1,i)= CCold(1,1,i)
1895 !c Ctilde(1,2,i)= CCold(1,2,i)
1896 !c Ctilde(2,1,i)=-CCold(2,1,i)
1897 !c Ctilde(2,2,i)=-CCold(2,2,i)
1899 !c Ctilde(1,1,i)=0.0d0
1900 !c Ctilde(1,2,i)=0.0d0
1901 !c Ctilde(2,1,i)=0.0d0
1902 !c Ctilde(2,2,i)=0.0d0
1903 DDold(1,1,i)= b(6,i)
1904 DDold(2,2,i)=-b(6,i)
1905 DDold(2,1,i)= b(8,i)
1906 DDold(1,2,i)= b(8,i)
1907 DDold(1,1,-i)= b(6,i)
1908 DDold(2,2,-i)=-b(6,i)
1909 DDold(2,1,-i)=-b(8,i)
1910 DDold(1,2,-i)=-b(8,i)
1915 !c Dtilde(1,1,i)= DD(1,1,i)
1916 !c Dtilde(1,2,i)= DD(1,2,i)
1917 !c Dtilde(2,1,i)=-DD(2,1,i)
1918 !c Dtilde(2,2,i)=-DD(2,2,i)
1920 !c Dtilde(1,1,i)=0.0d0
1921 !c Dtilde(1,2,i)=0.0d0
1922 !c Dtilde(2,1,i)=0.0d0
1923 !c Dtilde(2,2,i)=0.0d0
1924 EEold(1,1,i)= b(10,i)+b(11,i)
1925 EEold(2,2,i)=-b(10,i)+b(11,i)
1926 EEold(2,1,i)= b(12,i)-b(13,i)
1927 EEold(1,2,i)= b(12,i)+b(13,i)
1928 EEold(1,1,-i)= b(10,i)+b(11,i)
1929 EEold(2,2,-i)=-b(10,i)+b(11,i)
1930 EEold(2,1,-i)=-b(12,i)+b(13,i)
1931 EEold(1,2,-i)=-b(12,i)-b(13,i)
1932 write(iout,*) "TU DOCHODZE"
1938 !c ee(2,1,i)=ee(1,2,i)
1943 "Coefficients of the cumulants (independent of valence angles)"
1944 do i=-nloctyp+1,nloctyp-1
1945 write (iout,*) 'Type ',onelet(iloctyp(i))
1947 write(iout,'(2f10.5)') B(3,i),B(5,i)
1949 write(iout,'(2f10.5)') B(2,i),B(4,i)
1952 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1956 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1960 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1969 ! Read torsional parameters in old format
1971 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1973 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1974 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1975 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1977 !el from energy module--------
1978 allocate(v1(nterm_old,ntortyp,ntortyp))
1979 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1980 !el---------------------------
1985 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1991 write (iout,'(/a/)') 'Torsional constants:'
1994 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1995 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2001 ! Read torsional parameters
2003 IF (TOR_MODE.eq.0) THEN
2004 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2005 read (itorp,*,end=113,err=113) ntortyp
2006 !el from energy module---------
2007 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2008 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2010 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2011 allocate(vlor2(maxlor,ntortyp,ntortyp))
2012 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2013 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2015 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2016 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2017 !el---------------------------
2020 !el---------------------------
2022 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2024 itortyp(i)=-itortyp(-i)
2026 itortyp(ntyp1)=ntortyp
2027 itortyp(-ntyp1)=-ntortyp
2029 write (iout,*) 'ntortyp',ntortyp
2031 do j=-ntortyp+1,ntortyp-1
2032 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2034 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2035 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2038 do k=1,nterm(i,j,iblock)
2039 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2041 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2042 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2043 v0ij=v0ij+si*v1(k,i,j,iblock)
2045 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2046 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2047 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2049 do k=1,nlor(i,j,iblock)
2050 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2051 vlor2(k,i,j),vlor3(k,i,j)
2052 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2055 v0(-i,-j,iblock)=v0ij
2061 write (iout,'(/a/)') 'Torsional constants:'
2063 do i=-ntortyp,ntortyp
2064 do j=-ntortyp,ntortyp
2065 write (iout,*) 'ityp',i,' jtyp',j
2066 write (iout,*) 'Fourier constants'
2067 do k=1,nterm(i,j,iblock)
2068 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2071 write (iout,*) 'Lorenz constants'
2072 do k=1,nlor(i,j,iblock)
2073 write (iout,'(3(1pe15.5))') &
2074 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2080 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2082 ! 6/23/01 Read parameters for double torsionals
2084 !el from energy module------------
2085 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2086 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2087 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2088 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2089 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2090 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2091 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2092 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2093 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2094 !---------------------------------
2098 do j=-ntortyp+1,ntortyp-1
2099 do k=-ntortyp+1,ntortyp-1
2100 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2101 ! write (iout,*) "OK onelett",
2104 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2105 .or. t3.ne.toronelet(k)) then
2106 write (iout,*) "Error in double torsional parameter file",&
2109 call MPI_Finalize(Ierror)
2111 stop "Error in double torsional parameter file"
2113 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2114 ntermd_2(i,j,k,iblock)
2115 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2116 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2117 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2118 ntermd_1(i,j,k,iblock))
2119 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2120 ntermd_1(i,j,k,iblock))
2121 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2122 ntermd_1(i,j,k,iblock))
2123 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2124 ntermd_1(i,j,k,iblock))
2125 ! Martix of D parameters for one dimesional foureir series
2126 do l=1,ntermd_1(i,j,k,iblock)
2127 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2128 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2129 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2130 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2131 ! write(iout,*) "whcodze" ,
2132 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2134 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2135 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2136 v2s(m,l,i,j,k,iblock),&
2137 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2138 ! Martix of D parameters for two dimesional fourier series
2139 do l=1,ntermd_2(i,j,k,iblock)
2141 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2142 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2143 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2144 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2153 write (iout,*) 'Constants for double torsionals'
2156 do j=-ntortyp+1,ntortyp-1
2157 do k=-ntortyp+1,ntortyp-1
2158 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2159 ' nsingle',ntermd_1(i,j,k,iblock),&
2160 ' ndouble',ntermd_2(i,j,k,iblock)
2162 write (iout,*) 'Single angles:'
2163 do l=1,ntermd_1(i,j,k,iblock)
2164 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2165 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2166 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2167 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2170 write (iout,*) 'Pairs of angles:'
2171 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2172 do l=1,ntermd_2(i,j,k,iblock)
2173 write (iout,'(i5,20f10.5)') &
2174 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2177 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2178 do l=1,ntermd_2(i,j,k,iblock)
2179 write (iout,'(i5,20f10.5)') &
2180 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2181 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2189 ELSE IF (TOR_MODE.eq.1) THEN
2191 !C read valence-torsional parameters
2192 read (itorp,*,end=121,err=121) ntortyp
2194 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2196 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2197 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2200 itype2loc(i)=itortyp(i)
2204 itortyp(i)=-itortyp(-i)
2206 do i=-ntortyp+1,ntortyp-1
2207 do j=-ntortyp+1,ntortyp-1
2208 !C first we read the cos and sin gamma parameters
2209 read (itorp,'(13x,a)',end=121,err=121) string
2210 write (iout,*) i,j,string
2211 read (itorp,*,end=121,err=121) &
2212 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2213 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2214 do k=1,nterm_kcc(j,i)
2215 do l=1,nterm_kcc_Tb(j,i)
2216 do ll=1,nterm_kcc_Tb(j,i)
2217 read (itorp,*,end=121,err=121) ii,jj,kk, &
2218 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2226 !c AL 4/8/16: Calculate coefficients from one-body parameters
2228 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2229 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2230 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2231 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2232 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2235 print *,i,itortyp(i)
2236 itortyp(i)=itype2loc(i)
2239 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2241 do i=-ntortyp+1,ntortyp-1
2242 do j=-ntortyp+1,ntortyp-1
2245 do k=1,nterm_kcc_Tb(j,i)
2246 do l=1,nterm_kcc_Tb(j,i)
2247 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2248 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2249 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2250 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2253 do k=1,nterm_kcc_Tb(j,i)
2254 do l=1,nterm_kcc_Tb(j,i)
2256 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2257 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2258 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2259 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2261 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2262 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2263 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2264 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2268 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2272 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2277 if (tor_mode.gt.0 .and. lprint) then
2278 !c Print valence-torsional parameters
2279 write (iout,'(a)') &
2280 "Parameters of the valence-torsional potentials"
2281 do i=-ntortyp+1,ntortyp-1
2282 do j=-ntortyp+1,ntortyp-1
2283 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2284 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2285 do k=1,nterm_kcc(j,i)
2286 do l=1,nterm_kcc_Tb(j,i)
2287 do ll=1,nterm_kcc_Tb(j,i)
2288 write (iout,'(3i5,2f15.4)')&
2289 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2297 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2298 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2299 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2300 !el from energy module---------
2301 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2302 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2304 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2305 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2306 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2307 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2309 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2310 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2311 !el---------------------------
2314 !el--------------------
2315 read (itorp_nucl,*,end=113,err=113) &
2316 (itortyp_nucl(i),i=1,ntyp_molec(2))
2317 ! print *,itortyp_nucl(:)
2318 !c write (iout,*) 'ntortyp',ntortyp
2321 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2322 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2325 do k=1,nterm_nucl(i,j)
2326 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2327 v0ij=v0ij+si*v1_nucl(k,i,j)
2330 do k=1,nlor_nucl(i,j)
2331 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2332 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2333 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2339 ! Read of Side-chain backbone correlation parameters
2340 ! Modified 11 May 2012 by Adasko
2343 read (isccor,*,end=119,err=119) nsccortyp
2345 !el from module energy-------------
2346 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2347 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2348 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2349 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2350 !-----------------------------------
2352 !el from module energy-------------
2353 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2355 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2357 isccortyp(i)=-isccortyp(-i)
2359 iscprol=isccortyp(20)
2360 ! write (iout,*) 'ntortyp',ntortyp
2362 !c maxinter is maximum interaction sites
2363 !el from module energy---------
2364 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2365 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2366 -nsccortyp:nsccortyp))
2367 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2368 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2369 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2370 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2371 !-----------------------------------
2375 read (isccor,*,end=119,err=119) &
2376 nterm_sccor(i,j),nlor_sccor(i,j)
2382 nterm_sccor(-i,j)=nterm_sccor(i,j)
2383 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2384 nterm_sccor(i,-j)=nterm_sccor(i,j)
2385 do k=1,nterm_sccor(i,j)
2386 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2388 if (j.eq.iscprol) then
2389 if (i.eq.isccortyp(10)) then
2390 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2391 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2393 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2394 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2395 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2396 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2397 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2398 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2399 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2400 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2403 if (i.eq.isccortyp(10)) then
2404 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2405 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2407 if (j.eq.isccortyp(10)) then
2408 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2409 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2411 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2412 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2413 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2414 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2415 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2416 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2420 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2421 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2422 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2423 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2426 do k=1,nlor_sccor(i,j)
2427 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2428 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2429 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2430 (1+vlor3sccor(k,i,j)**2)
2432 v0sccor(l,i,j)=v0ijsccor
2433 v0sccor(l,-i,j)=v0ijsccor1
2434 v0sccor(l,i,-j)=v0ijsccor2
2435 v0sccor(l,-i,-j)=v0ijsccor3
2441 !el from module energy-------------
2442 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2444 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2445 ! write (iout,*) 'ntortyp',ntortyp
2447 !c maxinter is maximum interaction sites
2448 !el from module energy---------
2449 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2450 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2451 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2452 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2453 !-----------------------------------
2457 read (isccor,*,end=119,err=119) &
2458 nterm_sccor(i,j),nlor_sccor(i,j)
2462 do k=1,nterm_sccor(i,j)
2463 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2465 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2468 do k=1,nlor_sccor(i,j)
2469 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2470 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2471 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2472 (1+vlor3sccor(k,i,j)**2)
2474 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2483 write (iout,'(/a/)') 'Torsional constants:'
2486 write (iout,*) 'ityp',i,' jtyp',j
2487 write (iout,*) 'Fourier constants'
2488 do k=1,nterm_sccor(i,j)
2489 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2491 write (iout,*) 'Lorenz constants'
2492 do k=1,nlor_sccor(i,j)
2493 write (iout,'(3(1pe15.5))') &
2494 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2501 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2502 ! interaction energy of the Gly, Ala, and Pro prototypes.
2505 ! Read electrostatic-interaction parameters
2510 write (iout,'(/a)') 'Electrostatic interaction constants:'
2511 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2512 'IT','JT','APP','BPP','AEL6','AEL3'
2514 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2515 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2516 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2517 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2522 app (i,j)=epp(i,j)*rri*rri
2523 bpp (i,j)=-2.0D0*epp(i,j)*rri
2524 ael6(i,j)=elpp6(i,j)*4.2D0**6
2525 ael3(i,j)=elpp3(i,j)*4.2D0**3
2527 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2533 ! Read side-chain interaction parameters.
2535 !el from module energy - COMMON.INTERACT-------
2536 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2537 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2538 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2540 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2541 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2542 allocate(epslip(ntyp,ntyp))
2550 !--------------------------------
2552 read (isidep,*,end=117,err=117) ipot,expon
2553 if (ipot.lt.1 .or. ipot.gt.5) then
2554 write (iout,'(2a)') 'Error while reading SC interaction',&
2555 'potential file - unknown potential type.'
2557 call MPI_Finalize(Ierror)
2563 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2564 ', exponents are ',expon,2*expon
2565 ! goto (10,20,30,30,40) ipot
2567 !----------------------- LJ potential ---------------------------------
2569 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2570 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2571 (sigma0(i),i=1,ntyp)
2573 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2574 write (iout,'(a/)') 'The epsilon array:'
2575 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2576 write (iout,'(/a)') 'One-body parameters:'
2577 write (iout,'(a,4x,a)') 'residue','sigma'
2578 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2581 !----------------------- LJK potential --------------------------------
2583 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2584 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2585 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2587 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2588 write (iout,'(a/)') 'The epsilon array:'
2589 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2590 write (iout,'(/a)') 'One-body parameters:'
2591 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2592 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2596 !---------------------- GB or BP potential -----------------------------
2599 ! print *,"I AM in SCELE",scelemode
2600 if (scelemode.eq.0) then
2602 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2604 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2605 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2606 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2607 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2609 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2612 ! For the GB potential convert sigma'**2 into chi'
2615 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2619 write (iout,'(/a/)') 'Parameters of the BP potential:'
2620 write (iout,'(a/)') 'The epsilon array:'
2621 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2622 write (iout,'(/a)') 'One-body parameters:'
2623 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2625 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2626 chip(i),alp(i),i=1,ntyp)
2629 ! print *,ntyp,"NTYP"
2630 allocate(icharge(ntyp1))
2631 ! print *,ntyp,icharge(i)
2633 read (isidep,*) (icharge(i),i=1,ntyp)
2634 print *,ntyp,icharge(i)
2635 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2636 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2637 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2638 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2639 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2640 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2641 allocate(epsintab(ntyp,ntyp))
2642 allocate(dtail(2,ntyp,ntyp))
2643 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2644 allocate(wqdip(2,ntyp,ntyp))
2645 allocate(wstate(4,ntyp,ntyp))
2646 allocate(dhead(2,2,ntyp,ntyp))
2647 allocate(nstate(ntyp,ntyp))
2648 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2649 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2652 ! write (*,*) "Im in ALAB", i, " ", j
2654 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2655 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2656 chis(i,j),chis(j,i), &
2657 nstate(i,j),(wstate(k,i,j),k=1,4), &
2658 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2659 dtail(1,i,j),dtail(2,i,j), &
2660 epshead(i,j),sig0head(i,j), &
2661 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2662 alphapol(i,j),alphapol(j,i), &
2663 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
2664 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2670 sigma(i,j) = sigma(j,i)
2671 sigmap1(i,j)=sigmap1(j,i)
2672 sigmap2(i,j)=sigmap2(j,i)
2673 sigiso1(i,j)=sigiso1(j,i)
2674 sigiso2(i,j)=sigiso2(j,i)
2675 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2676 nstate(i,j) = nstate(j,i)
2677 dtail(1,i,j) = dtail(1,j,i)
2678 dtail(2,i,j) = dtail(2,j,i)
2680 alphasur(k,i,j) = alphasur(k,j,i)
2681 wstate(k,i,j) = wstate(k,j,i)
2682 alphiso(k,i,j) = alphiso(k,j,i)
2685 dhead(2,1,i,j) = dhead(1,1,j,i)
2686 dhead(2,2,i,j) = dhead(1,2,j,i)
2687 dhead(1,1,i,j) = dhead(2,1,j,i)
2688 dhead(1,2,i,j) = dhead(2,2,j,i)
2690 epshead(i,j) = epshead(j,i)
2691 sig0head(i,j) = sig0head(j,i)
2694 wqdip(k,i,j) = wqdip(k,j,i)
2697 wquad(i,j) = wquad(j,i)
2698 epsintab(i,j) = epsintab(j,i)
2699 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2704 !--------------------- GBV potential -----------------------------------
2706 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2707 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2708 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2709 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2711 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2712 write (iout,'(a/)') 'The epsilon array:'
2713 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2714 write (iout,'(/a)') 'One-body parameters:'
2715 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2716 's||/s_|_^2',' chip ',' alph '
2717 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2718 sigii(i),chip(i),alp(i),i=1,ntyp)
2721 write(iout,*)"Wrong ipot"
2727 !-----------------------------------------------------------------------
2728 ! Calculate the "working" parameters of SC interactions.
2730 !el from module energy - COMMON.INTERACT-------
2731 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2732 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2733 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2734 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2735 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2736 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2738 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2744 if (scelemode.eq.0) then
2753 sc_aa_tube_par(:)=0.0d0
2754 sc_bb_tube_par(:)=0.0d0
2756 !--------------------------------
2761 epslip(i,j)=epslip(j,i)
2764 if (scelemode.eq.0) then
2767 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2768 sigma(j,i)=sigma(i,j)
2769 rs0(i,j)=dwa16*sigma(i,j)
2774 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2775 'Working parameters of the SC interactions:',&
2776 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2781 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2783 ! print *,"SIGMA ZLA?",sigma(i,j)
2791 sigeps=dsign(1.0D0,epsij)
2793 aa_aq(i,j)=epsij*rrij*rrij
2794 ! print *,"ADASKO",epsij,rrij,expon
2795 bb_aq(i,j)=-sigeps*epsij*rrij
2796 aa_aq(j,i)=aa_aq(i,j)
2797 bb_aq(j,i)=bb_aq(i,j)
2798 epsijlip=epslip(i,j)
2799 sigeps=dsign(1.0D0,epsijlip)
2800 epsijlip=dabs(epsijlip)
2801 aa_lip(i,j)=epsijlip*rrij*rrij
2802 bb_lip(i,j)=-sigeps*epsijlip*rrij
2803 aa_lip(j,i)=aa_lip(i,j)
2804 bb_lip(j,i)=bb_lip(i,j)
2805 !C write(iout,*) aa_lip
2806 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2807 sigt1sq=sigma0(i)**2
2808 sigt2sq=sigma0(j)**2
2811 ratsig1=sigt2sq/sigt1sq
2812 ratsig2=1.0D0/ratsig1
2813 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2814 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2815 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2819 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2820 sigmaii(i,j)=rsum_max
2821 sigmaii(j,i)=rsum_max
2823 ! sigmaii(i,j)=r0(i,j)
2824 ! sigmaii(j,i)=r0(i,j)
2826 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2827 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2828 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2829 augm(i,j)=epsij*r_augm**(2*expon)
2830 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2837 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2838 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2839 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2844 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2845 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2846 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2847 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2848 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2849 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2850 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2851 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2852 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2853 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2854 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2855 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2856 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2857 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2858 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2859 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2868 read (isidep_nucl,*) ipot_nucl
2869 ! print *,"TU?!",ipot_nucl
2870 if (ipot_nucl.eq.1) then
2871 do i=1,ntyp_molec(2)
2872 do j=i,ntyp_molec(2)
2873 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2874 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2878 do i=1,ntyp_molec(2)
2879 do j=i,ntyp_molec(2)
2880 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2881 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2882 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2886 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2887 do i=1,ntyp_molec(2)
2888 do j=i,ntyp_molec(2)
2889 rrij=sigma_nucl(i,j)
2893 epsij=4*eps_nucl(i,j)
2894 sigeps=dsign(1.0D0,epsij)
2896 aa_nucl(i,j)=epsij*rrij*rrij
2897 bb_nucl(i,j)=-sigeps*epsij*rrij
2898 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2899 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2900 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2901 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2902 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2903 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2906 aa_nucl(i,j)=aa_nucl(j,i)
2907 bb_nucl(i,j)=bb_nucl(j,i)
2908 ael3_nucl(i,j)=ael3_nucl(j,i)
2909 ael6_nucl(i,j)=ael6_nucl(j,i)
2910 ael63_nucl(i,j)=ael63_nucl(j,i)
2911 ael32_nucl(i,j)=ael32_nucl(j,i)
2912 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2913 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2914 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2915 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2916 eps_nucl(i,j)=eps_nucl(j,i)
2917 sigma_nucl(i,j)=sigma_nucl(j,i)
2918 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2922 write(iout,*) "tube param"
2923 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2924 ccavtubpep,dcavtubpep,tubetranenepep
2925 sigmapeptube=sigmapeptube**6
2926 sigeps=dsign(1.0D0,epspeptube)
2927 epspeptube=dabs(epspeptube)
2928 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2929 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2930 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2932 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2933 ccavtub(i),dcavtub(i),tubetranene(i)
2934 sigmasctube=sigmasctube**6
2935 sigeps=dsign(1.0D0,epssctube)
2936 epssctube=dabs(epssctube)
2937 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2938 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2939 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2941 !-----------------READING SC BASE POTENTIALS-----------------------------
2942 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2943 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2944 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2945 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2946 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2947 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2948 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2949 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2950 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2951 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2952 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2953 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2954 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2955 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2956 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2957 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2960 do i=1,ntyp_molec(1)
2961 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2962 write (*,*) "Im in ", i, " ", j
2963 read(isidep_scbase,*) &
2964 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2965 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2966 write(*,*) "eps",eps_scbase(i,j)
2967 read(isidep_scbase,*) &
2968 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2969 chis_scbase(i,j,1),chis_scbase(i,j,2)
2970 read(isidep_scbase,*) &
2971 dhead_scbasei(i,j), &
2972 dhead_scbasej(i,j), &
2973 rborn_scbasei(i,j),rborn_scbasej(i,j)
2974 read(isidep_scbase,*) &
2975 (wdipdip_scbase(k,i,j),k=1,3), &
2976 (wqdip_scbase(k,i,j),k=1,2)
2977 read(isidep_scbase,*) &
2978 alphapol_scbase(i,j), &
2979 epsintab_scbase(i,j)
2982 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2983 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2985 do i=1,ntyp_molec(1)
2986 do j=1,ntyp_molec(2)-1
2987 epsij=eps_scbase(i,j)
2988 rrij=sigma_scbase(i,j)
2993 sigeps=dsign(1.0D0,epsij)
2995 aa_scbase(i,j)=epsij*rrij*rrij
2996 bb_scbase(i,j)=-sigeps*epsij*rrij
2999 !-----------------READING PEP BASE POTENTIALS-------------------
3000 allocate(eps_pepbase(ntyp_molec(2)))
3001 allocate(sigma_pepbase(ntyp_molec(2)))
3002 allocate(chi_pepbase(ntyp_molec(2),2))
3003 allocate(chipp_pepbase(ntyp_molec(2),2))
3004 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3005 allocate(sigmap1_pepbase(ntyp_molec(2)))
3006 allocate(sigmap2_pepbase(ntyp_molec(2)))
3007 allocate(chis_pepbase(ntyp_molec(2),2))
3008 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3011 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3012 write (*,*) "Im in ", i, " ", j
3013 read(isidep_pepbase,*) &
3014 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3015 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3016 write(*,*) "eps",eps_pepbase(j)
3017 read(isidep_pepbase,*) &
3018 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3019 chis_pepbase(j,1),chis_pepbase(j,2)
3020 read(isidep_pepbase,*) &
3021 (wdipdip_pepbase(k,j),k=1,3)
3023 allocate(aa_pepbase(ntyp_molec(2)))
3024 allocate(bb_pepbase(ntyp_molec(2)))
3026 do j=1,ntyp_molec(2)-1
3027 epsij=eps_pepbase(j)
3028 rrij=sigma_pepbase(j)
3033 sigeps=dsign(1.0D0,epsij)
3035 aa_pepbase(j)=epsij*rrij*rrij
3036 bb_pepbase(j)=-sigeps*epsij*rrij
3038 !--------------READING SC PHOSPHATE-------------------------------------
3039 allocate(eps_scpho(ntyp_molec(1)))
3040 allocate(sigma_scpho(ntyp_molec(1)))
3041 allocate(chi_scpho(ntyp_molec(1),2))
3042 allocate(chipp_scpho(ntyp_molec(1),2))
3043 allocate(alphasur_scpho(4,ntyp_molec(1)))
3044 allocate(sigmap1_scpho(ntyp_molec(1)))
3045 allocate(sigmap2_scpho(ntyp_molec(1)))
3046 allocate(chis_scpho(ntyp_molec(1),2))
3047 allocate(wqq_scpho(ntyp_molec(1)))
3048 allocate(wqdip_scpho(2,ntyp_molec(1)))
3049 allocate(alphapol_scpho(ntyp_molec(1)))
3050 allocate(epsintab_scpho(ntyp_molec(1)))
3051 allocate(dhead_scphoi(ntyp_molec(1)))
3052 allocate(rborn_scphoi(ntyp_molec(1)))
3053 allocate(rborn_scphoj(ntyp_molec(1)))
3054 allocate(alphi_scpho(ntyp_molec(1)))
3058 do j=1,ntyp_molec(1) ! without U then we will take T for U
3059 write (*,*) "Im in scpho ", i, " ", j
3060 read(isidep_scpho,*) &
3061 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3062 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3063 write(*,*) "eps",eps_scpho(j)
3064 read(isidep_scpho,*) &
3065 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3066 chis_scpho(j,1),chis_scpho(j,2)
3067 read(isidep_scpho,*) &
3068 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3069 read(isidep_scpho,*) &
3070 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3074 allocate(aa_scpho(ntyp_molec(1)))
3075 allocate(bb_scpho(ntyp_molec(1)))
3077 do j=1,ntyp_molec(1)
3084 sigeps=dsign(1.0D0,epsij)
3086 aa_scpho(j)=epsij*rrij*rrij
3087 bb_scpho(j)=-sigeps*epsij*rrij
3091 read(isidep_peppho,*) &
3092 eps_peppho,sigma_peppho
3093 read(isidep_peppho,*) &
3094 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3095 read(isidep_peppho,*) &
3096 (wqdip_peppho(k),k=1,2)
3104 sigeps=dsign(1.0D0,epsij)
3106 aa_peppho=epsij*rrij*rrij
3107 bb_peppho=-sigeps*epsij*rrij
3110 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3115 ! Define the SC-p interaction constants (hard-coded; old style)
3118 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3120 ! aad(i,1)=0.3D0*4.0D0**12
3121 ! Following line for constants currently implemented
3122 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3123 aad(i,1)=1.5D0*4.0D0**12
3124 ! aad(i,1)=0.17D0*5.6D0**12
3126 ! "Soft" SC-p repulsion
3128 ! Following line for constants currently implemented
3129 ! aad(i,1)=0.3D0*4.0D0**6
3130 ! "Hard" SC-p repulsion
3131 bad(i,1)=3.0D0*4.0D0**6
3132 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3141 ! 8/9/01 Read the SC-p interaction constants from file
3144 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3147 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3148 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3149 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3150 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3154 write (iout,*) "Parameters of SC-p interactions:"
3156 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3157 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3162 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3164 do i=1,ntyp_molec(2)
3165 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3167 do i=1,ntyp_molec(2)
3168 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3169 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3171 r0pp=1.12246204830937298142*5.16158
3177 ! Define the constants of the disulfide bridge
3181 ! Old arbitrary potential - commented out.
3186 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3187 ! energy surface of diethyl disulfide.
3188 ! A. Liwo and U. Kozlowska, 11/24/03
3205 write (iout,'(/a)') "Disulfide bridge parameters:"
3206 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3207 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3208 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3209 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3212 if (shield_mode.gt.0) then
3213 pi=4.0D0*datan(1.0D0)
3214 !C VSolvSphere the volume of solving sphere
3216 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3217 !C there will be no distinction between proline peptide group and normal peptide
3218 !C group in case of shielding parameters
3219 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3220 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3221 write (iout,*) VSolvSphere,VSolvSphere_div
3222 !C long axis of side chain
3224 long_r_sidechain(i)=vbldsc0(1,i)
3225 ! if (scelemode.eq.0) then
3226 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3227 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3229 ! short_r_sidechain(i)=sigma(i,i)
3231 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3238 111 write (iout,*) "Error reading bending energy parameters."
3240 112 write (iout,*) "Error reading rotamer energy parameters."
3242 113 write (iout,*) "Error reading torsional energy parameters."
3244 114 write (iout,*) "Error reading double torsional energy parameters."
3246 115 write (iout,*) &
3247 "Error reading cumulant (multibody energy) parameters."
3249 116 write (iout,*) "Error reading electrostatic energy parameters."
3251 117 write (iout,*) "Error reading side chain interaction parameters."
3253 118 write (iout,*) "Error reading SCp interaction parameters."
3255 119 write (iout,*) "Error reading SCCOR parameters"
3257 121 write (iout,*) "Error in Czybyshev parameters"
3260 call MPI_Finalize(Ierror)
3264 end subroutine parmread
3266 !-----------------------------------------------------------------------------
3268 !-----------------------------------------------------------------------------
3269 subroutine printmat(ldim,m,n,iout,key,a)
3272 character(len=3),dimension(n) :: key
3273 real(kind=8),dimension(ldim,n) :: a
3275 integer :: i,j,k,m,iout,nlim
3279 write (iout,1000) (key(k),k=i,nlim)
3281 1000 format (/5x,8(6x,a3))
3282 1020 format (/80(1h-)/)
3284 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3287 1010 format (a3,2x,8(f9.4))
3289 end subroutine printmat
3290 !-----------------------------------------------------------------------------
3292 !-----------------------------------------------------------------------------
3294 ! Read the PDB file and convert the peptide geometry into virtual-chain
3297 use energy_data, only: itype,istype
3301 ! use control, only: rescode,sugarcode
3302 ! implicit real*8 (a-h,o-z)
3303 ! include 'DIMENSIONS'
3304 ! include 'COMMON.LOCAL'
3305 ! include 'COMMON.VAR'
3306 ! include 'COMMON.CHAIN'
3307 ! include 'COMMON.INTERACT'
3308 ! include 'COMMON.IOUNITS'
3309 ! include 'COMMON.GEO'
3310 ! include 'COMMON.NAMES'
3311 ! include 'COMMON.CONTROL'
3312 ! include 'COMMON.DISTFIT'
3313 ! include 'COMMON.SETUP'
3314 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3316 logical :: lprn=.true.,fail
3317 real(kind=8),dimension(3) :: e1,e2,e3
3318 real(kind=8) :: dcj,efree_temp
3319 character(len=3) :: seq,res,res2
3320 character(len=5) :: atom
3321 character(len=80) :: card
3322 real(kind=8),dimension(3,20) :: sccor
3323 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3324 integer :: isugar,molecprev,firstion
3325 character*1 :: sugar
3327 real(kind=8),dimension(3) :: ccc
3329 integer,dimension(2,maxres/3) :: hfrag_alloc
3330 integer,dimension(4,maxres/3) :: bfrag_alloc
3331 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3332 real(kind=8),dimension(:,:), allocatable :: c_temporary
3333 integer,dimension(:,:) , allocatable :: itype_temporary
3334 integer,dimension(:),allocatable :: istype_temp
3341 ! write (2,*) "UNRES_PDB",unres_pdb
3361 !-----------------------------
3362 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3363 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3364 if(.not. allocated(istype)) allocate(istype(maxres))
3366 read (ipdbin,'(a80)',end=10) card
3367 write (iout,'(a)') card
3368 if (card(:5).eq.'HELIX') then
3371 read(card(22:25),*) hfrag(1,nhfrag)
3372 read(card(34:37),*) hfrag(2,nhfrag)
3374 if (card(:5).eq.'SHEET') then
3377 read(card(24:26),*) bfrag(1,nbfrag)
3378 read(card(35:37),*) bfrag(2,nbfrag)
3379 !rc----------------------------------------
3380 !rc to be corrected !!!
3381 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3382 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3383 !rc----------------------------------------
3385 if (card(:3).eq.'END') then
3387 else if (card(:3).eq.'TER') then
3392 itype(ires_old,molecule)=ntyp1_molec(molecule)
3393 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3394 nres_molec(molecule)=nres_molec(molecule)+2
3396 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3399 dc(j,ires)=sccor(j,iii)
3402 call sccenter(ires,iii,sccor)
3408 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3409 ! Fish out the ATOM cards.
3410 ! write(iout,*) 'card',card(1:20)
3411 ! print *,"ATU ",card(1:6), CARD(3:6)
3413 if (index(card(1:4),'ATOM').gt.0) then
3414 read (card(12:16),*) atom
3415 ! write (iout,*) "! ",atom," !",ires
3416 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3417 read (card(23:26),*) ires
3418 read (card(18:20),'(a3)') res
3419 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3420 ! & " ires_old",ires_old
3421 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3422 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3423 if (ires-ishift+ishift1.ne.ires_old) then
3424 ! Calculate the CM of the preceding residue.
3425 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3427 ! write (iout,*) "Calculating sidechain center iii",iii
3430 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3433 call sccenter(ires_old,iii,sccor)
3437 ! Start new residue.
3438 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3441 else if (ibeg.eq.1) then
3442 write (iout,*) "BEG ires",ires
3444 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3447 nres_molec(molecule)=nres_molec(molecule)+1
3449 ires=ires-ishift+ishift1
3451 ! write (iout,*) "ishift",ishift," ires",ires,&
3452 ! " ires_old",ires_old
3454 else if (ibeg.eq.2) then
3456 ishift=-ires_old+ires-1 !!!!!
3457 ishift1=ishift1-1 !!!!!
3458 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3459 ires=ires-ishift+ishift1
3460 ! print *,ires,ishift,ishift1
3464 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3465 ires=ires-ishift+ishift1
3468 ! print *,'atom',ires,atom
3469 if (res.eq.'ACE' .or. res.eq.'NHE') then
3472 if (atom.eq.'CA '.or.atom.eq.'N ') then
3474 itype(ires,molecule)=rescode(ires,res,0,molecule)
3476 ! nres_molec(molecule)=nres_molec(molecule)+1
3480 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3481 ! nres_molec(molecule)=nres_molec(molecule)+1
3482 read (card(19:19),'(a1)') sugar
3483 isugar=sugarcode(sugar,ires)
3484 ! if (ibeg.eq.1) then
3488 ! print *,"ires=",ires,istype(ires)
3494 ires=ires-ishift+ishift1
3496 ! write (iout,*) "ires_old",ires_old," ires",ires
3497 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3500 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3501 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3502 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3503 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3504 ! print *,ires,ishift,ishift1
3505 ! write (iout,*) "backbone ",atom
3507 write (iout,'(2i3,2x,a,3f8.3)') &
3508 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3511 nres_molec(molecule)=nres_molec(molecule)+1
3513 sccor(j,iii)=c(j,ires)
3515 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3516 atom.eq."C2'" .or. atom.eq."C3'" &
3517 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3518 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3519 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3520 ! print *,ires,ishift,ishift1
3524 ! sccor(j,iii)=c(j,ires)
3527 c(j,ires)=c(j,ires)+ccc(j)/5.0
3529 print *,counter,molecule
3530 if (counter.eq.5) then
3532 nres_molec(molecule)=nres_molec(molecule)+1
3535 ! sccor(j,iii)=c(j,ires)
3539 ! print *, "ATOM",atom(1:3)
3540 else if (atom.eq."C5'") then
3541 read (card(19:19),'(a1)') sugar
3542 isugar=sugarcode(sugar,ires)
3547 ! print *,ires,istype(ires)
3551 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3552 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3553 nres_molec(molecule)=nres_molec(molecule)+1
3554 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3558 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3560 else if ((atom.eq."C1'").and.unres_pdb) then
3562 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3563 ! write (*,*) card(23:27),ires,itype(ires,1)
3564 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3565 atom.ne.'N' .and. atom.ne.'C' .and. &
3566 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3567 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3568 .and. atom.ne.'P '.and. &
3569 atom(1:1).ne.'H' .and. &
3570 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3572 ! write (iout,*) "sidechain ",atom
3573 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3574 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3575 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3577 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3580 ! print *,"IONS",ions,card(1:6)
3581 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3582 if (firstion.eq.0) then
3586 dc(j,ires)=sccor(j,iii)
3589 call sccenter(ires,iii,sccor)
3592 read (card(12:16),*) atom
3593 ! print *,"HETATOM", atom
3594 read (card(18:20),'(a3)') res
3595 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3596 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3597 .or.(atom(1:2).eq.'K ')) &
3600 if (molecule.ne.5) molecprev=molecule
3602 nres_molec(molecule)=nres_molec(molecule)+1
3603 print *,"HERE",nres_molec(molecule)
3605 itype(ires,molecule)=rescode(ires,res,0,molecule)
3606 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3610 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3611 if (ires.eq.0) return
3612 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3615 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3617 nres_molec(molecule)=nres_molec(molecule)-2
3618 print *,'I have',nres, nres_molec(:)
3620 do k=1,4 ! ions are without dummy
3621 if (nres_molec(k).eq.0) cycle
3623 ! write (iout,*) i,itype(i,1)
3624 ! if (itype(i,1).eq.ntyp1) then
3625 ! write (iout,*) "dummy",i,itype(i,1)
3627 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3628 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3632 if (itype(i,k).eq.ntyp1_molec(k)) then
3633 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3634 if (itype(i+2,k).eq.0) then
3635 ! print *,"masz sieczke"
3637 if (itype(i+2,j).ne.0) then
3639 itype(i+1,j)=ntyp1_molec(j)
3640 nres_molec(k)=nres_molec(k)-1
3641 nres_molec(j)=nres_molec(j)+1
3647 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3648 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3649 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3650 ! if (unres_pdb) then
3651 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3652 ! print *,i,'tu dochodze'
3653 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3661 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3665 dcj=(c(j,i-2)-c(j,i-3))/2.0
3666 if (dcj.eq.0) dcj=1.23591524223
3671 else !itype(i+1,1).eq.ntyp1
3672 ! if (unres_pdb) then
3673 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3674 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3681 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
3682 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3686 dcj=(c(j,i+3)-c(j,i+2))/2.0
3687 if (dcj.eq.0) dcj=1.23591524223
3692 endif !itype(i+1,1).eq.ntyp1
3693 endif !itype.eq.ntyp1
3697 ! Calculate the CM of the last side chain.
3701 dc(j,ires)=sccor(j,iii)
3704 call sccenter(ires,iii,sccor)
3710 ! print *,"molecule",molecule
3711 if ((itype(nres,1).ne.10)) then
3713 if (molecule.eq.5) molecule=molecprev
3714 itype(nres,molecule)=ntyp1_molec(molecule)
3715 nres_molec(molecule)=nres_molec(molecule)+1
3716 ! if (unres_pdb) then
3717 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3718 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3725 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3729 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3730 c(j,nres)=c(j,nres-1)+dcj
3731 c(j,2*nres)=c(j,nres)
3735 ! print *,'I have',nres, nres_molec(:)
3737 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3739 if (nres.ne.nres0) then
3740 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3742 stop "Error nres value in WHAM input"
3745 !---------------------------------
3746 !el reallocate tables
3749 ! hfrag_alloc(j,i)=hfrag(j,i)
3752 ! bfrag_alloc(j,i)=bfrag(j,i)
3758 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3759 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3760 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3761 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3765 ! hfrag(j,i)=hfrag_alloc(j,i)
3770 ! bfrag(j,i)=bfrag_alloc(j,i)
3773 !el end reallocate tables
3774 !---------------------------------
3782 c(j,2*nres)=c(j,nres)
3785 if (itype(1,1).eq.ntyp1) then
3789 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3790 call refsys(2,3,4,e1,e2,e3,fail)
3797 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3798 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3802 dcj=(c(j,4)-c(j,3))/2.0
3808 ! First lets assign correct dummy to correct type of chain
3810 if (itype(1,1).eq.ntyp1) then
3811 if (itype(2,1).eq.0) then
3813 if (itype(2,j).ne.0) then
3815 itype(1,j)=ntyp1_molec(j)
3816 nres_molec(1)=nres_molec(1)-1
3817 nres_molec(j)=nres_molec(j)+1
3824 print *,'I have',nres, nres_molec(:)
3826 ! Copy the coordinates to reference coordinates
3832 ! Calculate internal coordinates.
3834 write (iout,'(/a)') &
3835 "Cartesian coordinates of the reference structure"
3836 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3837 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3839 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3840 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3841 (c(j,ires+nres),j=1,3)
3844 ! znamy już nres wiec mozna alokowac tablice
3845 ! Calculate internal coordinates.
3846 if(me.eq.king.or..not.out1file)then
3847 write (iout,'(a)') &
3848 "Backbone and SC coordinates as read from the PDB"
3850 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
3851 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3852 (c(j,nres+ires),j=1,3)
3855 ! NOW LETS ROCK! SORTING
3856 allocate(c_temporary(3,2*nres))
3857 allocate(itype_temporary(nres,5))
3858 if (.not.allocated(molnum)) allocate(molnum(nres+1))
3859 if (.not.allocated(istype)) write(iout,*) &
3860 "SOMETHING WRONG WITH ISTYTPE"
3861 allocate(istype_temp(nres))
3862 itype_temporary(:,:)=0
3866 if (itype(i,k).ne.0) then
3868 c_temporary(j,seqalingbegin)=c(j,i)
3869 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3872 itype_temporary(seqalingbegin,k)=itype(i,k)
3873 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3874 istype_temp(seqalingbegin)=istype(i)
3875 molnum(seqalingbegin)=k
3876 seqalingbegin=seqalingbegin+1
3882 c(j,i)=c_temporary(j,i)
3887 itype(i,k)=itype_temporary(i,k)
3888 istype(i)=istype_temp(i)
3891 ! if (itype(1,1).eq.ntyp1) then
3894 ! if (unres_pdb) then
3895 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3896 ! call refsys(2,3,4,e1,e2,e3,fail)
3903 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3907 ! dcj=(c(j,4)-c(j,3))/2.0
3909 ! c(j,nres+1)=c(j,1)
3915 write (iout,'(/a)') &
3916 "Cartesian coordinates of the reference structure after sorting"
3917 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3918 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3920 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3921 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3922 (c(j,ires+nres),j=1,3)
3926 ! print *,seqalingbegin,nres
3927 if(.not.allocated(vbld)) then
3928 allocate(vbld(2*nres))
3933 if(.not.allocated(vbld_inv)) then
3934 allocate(vbld_inv(2*nres))
3940 if(.not.allocated(theta)) then
3941 allocate(theta(nres+2))
3945 if(.not.allocated(phi)) allocate(phi(nres+2))
3946 if(.not.allocated(alph)) allocate(alph(nres+2))
3947 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3948 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3949 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3950 if(.not.allocated(costtab)) allocate(costtab(nres))
3951 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3952 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3953 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3954 if(.not.allocated(xxref)) allocate(xxref(nres))
3955 if(.not.allocated(yyref)) allocate(yyref(nres))
3956 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3957 if(.not.allocated(dc_norm)) then
3958 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3959 allocate(dc_norm(3,0:2*nres+2))
3963 call int_from_cart(.true.,.false.)
3964 call sc_loc_geom(.false.)
3966 thetaref(i)=theta(i)
3976 dc(j,i)=c(j,i+1)-c(j,i)
3977 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3982 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3983 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3985 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3989 ! Copy the coordinates to reference coordinates
3990 ! Splits to single chain if occurs
3994 ! cref(j,i,cou)=c(j,i)
3998 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3999 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4000 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4001 !-----------------------------
4005 write (iout,*) "symetr", symetr
4008 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4010 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4013 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4018 cref(j,i,cou)=c(j,i)
4019 cref(j,i+nres,cou)=c(j,i+nres)
4021 chain_rep(j,lll,kkk)=c(j,i)
4022 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4026 write (iout,*) chain_length
4027 if (chain_length.eq.0) chain_length=nres
4029 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4030 chain_rep(j,chain_length+nres,symetr) &
4031 =chain_rep(j,chain_length+nres,1)
4034 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4036 ! do kkk=1,chain_length
4037 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4041 ! makes copy of chains
4042 write (iout,*) "symetr", symetr
4047 if (symetr.gt.1) then
4054 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4060 write (iout,*) i,icha
4061 do lll=1,chain_length
4063 if (cou.le.nres) then
4065 kupa=mod(lll,chain_length)
4066 iprzes=(kkk-1)*chain_length+lll
4067 if (kupa.eq.0) kupa=chain_length
4068 write (iout,*) "kupa", kupa
4069 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4070 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4077 !-koniec robienia kopii
4080 write (iout,*) "nowa struktura", nperm
4082 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4084 cref(3,i,kkk),cref(1,nres+i,kkk),&
4085 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4087 100 format (//' alpha-carbon coordinates ',&
4088 ' centroid coordinates'/ &
4089 ' ', 6X,'X',11X,'Y',11X,'Z', &
4090 10X,'X',11X,'Y',11X,'Z')
4091 110 format (a,'(',i5,')',6f12.5)
4097 bfrag(i,j)=bfrag(i,j)-ishift
4103 hfrag(i,j)=hfrag(i,j)-ishift
4109 end subroutine readpdb
4110 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4111 !-----------------------------------------------------------------------------
4113 !-----------------------------------------------------------------------------
4114 subroutine read_control
4128 use random, only: random_init
4129 ! implicit real*8 (a-h,o-z)
4130 ! include 'DIMENSIONS'
4132 use prng, only:prng_restart
4134 logical :: OKRandom!, prng_restart
4137 ! include 'COMMON.IOUNITS'
4138 ! include 'COMMON.TIME1'
4139 ! include 'COMMON.THREAD'
4140 ! include 'COMMON.SBRIDGE'
4141 ! include 'COMMON.CONTROL'
4142 ! include 'COMMON.MCM'
4143 ! include 'COMMON.MAP'
4144 ! include 'COMMON.HEADER'
4145 ! include 'COMMON.CSA'
4146 ! include 'COMMON.CHAIN'
4147 ! include 'COMMON.MUCA'
4148 ! include 'COMMON.MD'
4149 ! include 'COMMON.FFIELD'
4150 ! include 'COMMON.INTERACT'
4151 ! include 'COMMON.SETUP'
4152 !el integer :: KDIAG,ICORFL,IXDR
4153 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4154 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4155 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4156 ! character(len=80) :: ucase
4157 character(len=640) :: controlcard
4159 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4165 read (INP,'(a)') titel
4166 call card_concat(controlcard,.true.)
4167 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4168 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4169 call reada(controlcard,'SEED',seed,0.0D0)
4170 call random_init(seed)
4171 ! Set up the time limit (caution! The time must be input in minutes!)
4172 read_cart=index(controlcard,'READ_CART').gt.0
4173 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4174 call readi(controlcard,'SYM',symetr,1)
4175 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4176 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4177 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4178 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4179 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4180 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4181 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4182 call reada(controlcard,'DRMS',drms,0.1D0)
4183 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4184 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4185 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4186 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4187 write (iout,'(a,f10.1)')'DRMS = ',drms
4188 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4189 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4191 call readi(controlcard,'NZ_START',nz_start,0)
4192 call readi(controlcard,'NZ_END',nz_end,0)
4193 call readi(controlcard,'IZ_SC',iz_sc,0)
4194 timlim=60.0D0*timlim
4195 safety = 60.0d0*safety
4198 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4199 !C SHIELD keyword sets if the shielding effect of side-chains is used
4200 !C 0 denots no shielding is used all peptide are equally despite the
4201 !C solvent accesible area
4202 !C 1 the newly introduced function
4203 !C 2 reseved for further possible developement
4204 call readi(controlcard,'SHIELD',shield_mode,0)
4205 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4206 write(iout,*) "shield_mode",shield_mode
4207 call readi(controlcard,'TORMODE',tor_mode,0)
4208 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4209 write(iout,*) "torsional and valence angle mode",tor_mode
4211 !C Varibles set size of box
4212 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4213 protein=index(controlcard,"PROTEIN").gt.0
4214 ions=index(controlcard,"IONS").gt.0
4215 nucleic=index(controlcard,"NUCLEIC").gt.0
4216 write (iout,*) "with_theta_constr ",with_theta_constr
4217 AFMlog=(index(controlcard,'AFM'))
4218 selfguide=(index(controlcard,'SELFGUIDE'))
4219 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4220 call readi(controlcard,'GENCONSTR',genconstr,0)
4221 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4222 call reada(controlcard,'BOXY',boxysize,100.0d0)
4223 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4224 call readi(controlcard,'TUBEMOD',tubemode,0)
4225 print *,"SCELE",scelemode
4226 call readi(controlcard,"SCELEMODE",scelemode,0)
4227 print *,"SCELE",scelemode
4229 ! elemode = 0 is orignal UNRES electrostatics
4230 ! elemode = 1 is "Momo" potentials in progress
4231 ! elemode = 2 is in development EVALD
4232 write (iout,*) TUBEmode,"TUBEMODE"
4233 if (TUBEmode.gt.0) then
4234 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4235 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4236 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4237 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4238 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4239 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4240 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4241 buftubebot=bordtubebot+tubebufthick
4242 buftubetop=bordtubetop-tubebufthick
4245 ! CUTOFFF ON ELECTROSTATICS
4246 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
4247 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4248 write(iout,*) "R_CUT_ELE=",r_cut_ele
4249 ! Lipidic parameters
4250 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4251 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4252 if (lipthick.gt.0.0d0) then
4253 bordliptop=(boxzsize+lipthick)/2.0
4254 bordlipbot=bordliptop-lipthick
4255 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4256 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4257 buflipbot=bordlipbot+lipbufthick
4258 bufliptop=bordliptop-lipbufthick
4259 if ((lipbufthick*2.0d0).gt.lipthick) &
4260 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4261 endif !lipthick.gt.0
4262 write(iout,*) "bordliptop=",bordliptop
4263 write(iout,*) "bordlipbot=",bordlipbot
4264 write(iout,*) "bufliptop=",bufliptop
4265 write(iout,*) "buflipbot=",buflipbot
4266 write (iout,*) "SHIELD MODE",shield_mode
4268 !C-------------------------
4269 minim=(index(controlcard,'MINIMIZE').gt.0)
4270 dccart=(index(controlcard,'CART').gt.0)
4271 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4272 overlapsc=.not.overlapsc
4273 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4274 searchsc=.not.searchsc
4275 sideadd=(index(controlcard,'SIDEADD').gt.0)
4276 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4277 outpdb=(index(controlcard,'PDBOUT').gt.0)
4278 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4279 pdbref=(index(controlcard,'PDBREF').gt.0)
4280 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4281 indpdb=index(controlcard,'PDBSTART')
4282 extconf=(index(controlcard,'EXTCONF').gt.0)
4283 call readi(controlcard,'IPRINT',iprint,0)
4284 call readi(controlcard,'MAXGEN',maxgen,10000)
4285 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4286 call readi(controlcard,"KDIAG",kdiag,0)
4287 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4288 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4289 write (iout,*) "RESCALE_MODE",rescale_mode
4290 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4291 if (index(controlcard,'REGULAR').gt.0.0D0) then
4292 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4296 if (index(controlcard,'CHECKGRAD').gt.0) then
4298 if (index(controlcard,'CART').gt.0) then
4300 elseif (index(controlcard,'CARINT').gt.0) then
4305 elseif (index(controlcard,'THREAD').gt.0) then
4307 call readi(controlcard,'THREAD',nthread,0)
4308 if (nthread.gt.0) then
4309 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4312 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4313 stop 'Error termination in Read_Control.'
4315 else if (index(controlcard,'MCMA').gt.0) then
4317 else if (index(controlcard,'MCEE').gt.0) then
4319 else if (index(controlcard,'MULTCONF').gt.0) then
4321 else if (index(controlcard,'MAP').gt.0) then
4323 call readi(controlcard,'MAP',nmap,0)
4324 else if (index(controlcard,'CSA').gt.0) then
4326 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4328 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4331 !fcm else if (index(controlcard,'MCMF').gt.0) then
4333 else if (index(controlcard,'SOFTREG').gt.0) then
4335 else if (index(controlcard,'CHECK_BOND').gt.0) then
4337 else if (index(controlcard,'TEST').gt.0) then
4339 else if (index(controlcard,'MD').gt.0) then
4341 else if (index(controlcard,'RE ').gt.0) then
4345 lmuca=index(controlcard,'MUCA').gt.0
4346 call readi(controlcard,'MUCADYN',mucadyn,0)
4347 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4348 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4350 write (iout,*) 'MUCADYN=',mucadyn
4351 write (iout,*) 'MUCASMOOTH=',muca_smooth
4354 iscode=index(controlcard,'ONE_LETTER')
4355 indphi=index(controlcard,'PHI')
4356 indback=index(controlcard,'BACK')
4357 iranconf=index(controlcard,'RAND_CONF')
4358 i2ndstr=index(controlcard,'USE_SEC_PRED')
4359 gradout=index(controlcard,'GRADOUT').gt.0
4360 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4361 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4362 if (me.eq.king .or. .not.out1file ) &
4363 write (iout,*) "DISTCHAINMAX",distchainmax
4365 if(me.eq.king.or..not.out1file) &
4366 write (iout,'(2a)') diagmeth(kdiag),&
4367 ' routine used to diagonalize matrices.'
4368 if (shield_mode.gt.0) then
4369 pi=4.0D0*datan(1.0D0)
4370 !C VSolvSphere the volume of solving sphere
4372 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4373 !C there will be no distinction between proline peptide group and normal peptide
4374 !C group in case of shielding parameters
4375 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4376 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4377 write (iout,*) VSolvSphere,VSolvSphere_div
4378 !C long axis of side chain
4380 ! long_r_sidechain(i)=vbldsc0(1,i)
4381 ! short_r_sidechain(i)=sigma0(i)
4382 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4388 end subroutine read_control
4389 !-----------------------------------------------------------------------------
4390 subroutine read_REMDpar
4392 ! Read REMD settings
4399 use control_data, only:out1file
4400 ! implicit real*8 (a-h,o-z)
4401 ! include 'DIMENSIONS'
4402 ! include 'COMMON.IOUNITS'
4403 ! include 'COMMON.TIME1'
4404 ! include 'COMMON.MD'
4407 !el include 'COMMON.LANGEVIN'
4409 !el include 'COMMON.LANGEVIN.lang0'
4411 ! include 'COMMON.INTERACT'
4412 ! include 'COMMON.NAMES'
4413 ! include 'COMMON.GEO'
4414 ! include 'COMMON.REMD'
4415 ! include 'COMMON.CONTROL'
4416 ! include 'COMMON.SETUP'
4417 ! character(len=80) :: ucase
4418 character(len=320) :: controlcard
4419 character(len=3200) :: controlcard1
4420 integer :: iremd_m_total
4423 ! real(kind=8) :: var,ene
4425 if(me.eq.king.or..not.out1file) &
4426 write (iout,*) "REMD setup"
4428 call card_concat(controlcard,.true.)
4429 call readi(controlcard,"NREP",nrep,3)
4430 call readi(controlcard,"NSTEX",nstex,1000)
4431 call reada(controlcard,"RETMIN",retmin,10.0d0)
4432 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4433 mremdsync=(index(controlcard,'SYNC').gt.0)
4434 call readi(controlcard,"NSYN",i_sync_step,100)
4435 restart1file=(index(controlcard,'REST1FILE').gt.0)
4436 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4437 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4438 if(max_cache_traj_use.gt.max_cache_traj) &
4439 max_cache_traj_use=max_cache_traj
4440 if(me.eq.king.or..not.out1file) then
4441 !d if (traj1file) then
4442 !rc caching is in testing - NTWX is not ignored
4443 !d write (iout,*) "NTWX value is ignored"
4444 !d write (iout,*) " trajectory is stored to one file by master"
4445 !d write (iout,*) " before exchange at NSTEX intervals"
4447 write (iout,*) "NREP= ",nrep
4448 write (iout,*) "NSTEX= ",nstex
4449 write (iout,*) "SYNC= ",mremdsync
4450 write (iout,*) "NSYN= ",i_sync_step
4451 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4454 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4455 if (index(controlcard,'TLIST').gt.0) then
4457 call card_concat(controlcard1,.true.)
4458 read(controlcard1,*) (remd_t(i),i=1,nrep)
4459 if(me.eq.king.or..not.out1file) &
4460 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4463 if (index(controlcard,'MLIST').gt.0) then
4465 call card_concat(controlcard1,.true.)
4466 read(controlcard1,*) (remd_m(i),i=1,nrep)
4467 if(me.eq.king.or..not.out1file) then
4468 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4471 iremd_m_total=iremd_m_total+remd_m(i)
4473 write (iout,*) 'Total number of replicas ',iremd_m_total
4476 if(me.eq.king.or..not.out1file) &
4477 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4479 end subroutine read_REMDpar
4480 !-----------------------------------------------------------------------------
4481 subroutine read_MDpar
4485 use control_data, only: r_cut,rlamb,out1file
4487 use geometry_data, only: pi
4489 ! implicit real*8 (a-h,o-z)
4490 ! include 'DIMENSIONS'
4491 ! include 'COMMON.IOUNITS'
4492 ! include 'COMMON.TIME1'
4493 ! include 'COMMON.MD'
4496 !el include 'COMMON.LANGEVIN'
4498 !el include 'COMMON.LANGEVIN.lang0'
4500 ! include 'COMMON.INTERACT'
4501 ! include 'COMMON.NAMES'
4502 ! include 'COMMON.GEO'
4503 ! include 'COMMON.SETUP'
4504 ! include 'COMMON.CONTROL'
4505 ! include 'COMMON.SPLITELE'
4506 ! character(len=80) :: ucase
4507 character(len=320) :: controlcard
4512 call card_concat(controlcard,.true.)
4513 call readi(controlcard,"NSTEP",n_timestep,1000000)
4514 call readi(controlcard,"NTWE",ntwe,100)
4515 call readi(controlcard,"NTWX",ntwx,1000)
4516 call reada(controlcard,"DT",d_time,1.0d-1)
4517 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4518 call reada(controlcard,"DAMAX",damax,1.0d1)
4519 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4520 call readi(controlcard,"LANG",lang,0)
4521 RESPA = index(controlcard,"RESPA") .gt. 0
4522 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4523 ntime_split0=ntime_split
4524 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4525 ntime_split0=ntime_split
4526 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4527 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4528 rest = index(controlcard,"REST").gt.0
4529 tbf = index(controlcard,"TBF").gt.0
4530 usampl = index(controlcard,"USAMPL").gt.0
4531 mdpdb = index(controlcard,"MDPDB").gt.0
4532 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4533 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4534 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4535 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4536 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4537 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4538 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4539 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4540 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4541 large = index(controlcard,"LARGE").gt.0
4542 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4543 rattle = index(controlcard,"RATTLE").gt.0
4544 preminim=(index(controlcard,'PREMINIM').gt.0)
4545 write (iout,*) "PREMINIM ",preminim
4546 dccart=(index(controlcard,'CART').gt.0)
4547 if (preminim) call read_minim
4548 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4554 if(me.eq.king.or..not.out1file) then
4556 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4558 write (iout,'(a)') "The units are:"
4559 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4560 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4561 " acceleration: angstrom/(48.9 fs)**2"
4562 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4564 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4565 write (iout,'(a60,f10.5,a)') &
4566 "Initial time step of numerical integration:",d_time,&
4568 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4570 write (iout,'(2a,i4,a)') &
4571 "A-MTS algorithm used; initial time step for fast-varying",&
4572 " short-range forces split into",ntime_split," steps."
4573 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4574 r_cut," lambda",rlamb
4576 write (iout,'(2a,f10.5)') &
4577 "Maximum acceleration threshold to reduce the time step",&
4578 "/increase split number:",damax
4579 write (iout,'(2a,f10.5)') &
4580 "Maximum predicted energy drift to reduce the timestep",&
4581 "/increase split number:",edriftmax
4582 write (iout,'(a60,f10.5)') &
4583 "Maximum velocity threshold to reduce velocities:",dvmax
4584 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4585 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4586 if (rattle) write (iout,'(a60)') &
4587 "Rattle algorithm used to constrain the virtual bonds"
4591 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4592 call reada(controlcard,"RWAT",rwat,1.4d0)
4593 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4594 surfarea=index(controlcard,"SURFAREA").gt.0
4595 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4596 if(me.eq.king.or..not.out1file)then
4597 write (iout,'(/a,$)') "Langevin dynamics calculation"
4599 write (iout,'(a/)') &
4600 " with direct integration of Langevin equations"
4601 else if (lang.eq.2) then
4602 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4603 else if (lang.eq.3) then
4604 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4605 else if (lang.eq.4) then
4606 write (iout,'(a/)') " in overdamped mode"
4608 write (iout,'(//a,i5)') &
4609 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4612 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4613 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4614 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4615 write (iout,'(a60,f10.5)') &
4616 "Scaling factor of the friction forces:",scal_fric
4617 if (surfarea) write (iout,'(2a,i10,a)') &
4618 "Friction coefficients will be scaled by solvent-accessible",&
4619 " surface area every",reset_fricmat," steps."
4621 ! Calculate friction coefficients and bounds of stochastic forces
4622 eta=6*pi*cPoise*etawat
4623 if(me.eq.king.or..not.out1file) &
4624 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4627 do j=1,5 !types of molecules
4628 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4629 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4631 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4632 do j=1,5 !types of molecules
4634 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4635 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4639 if(me.eq.king.or..not.out1file)then
4640 write (iout,'(/2a/)') &
4641 "Radii of site types and friction coefficients and std's of",&
4642 " stochastic forces of fully exposed sites"
4643 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4645 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4646 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4650 if(me.eq.king.or..not.out1file)then
4651 write (iout,'(a)') "Berendsen bath calculation"
4652 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4653 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4655 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4656 count_reset_moment," steps"
4658 write (iout,'(a,i10,a)') &
4659 "Velocities will be reset at random every",count_reset_vel,&
4663 if(me.eq.king.or..not.out1file) &
4664 write (iout,'(a31)') "Microcanonical mode calculation"
4666 if(me.eq.king.or..not.out1file)then
4667 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4669 write(iout,*) "MD running with constraints."
4670 write(iout,*) "Equilibration time ", eq_time, " mtus."
4671 write(iout,*) "Constraining ", nfrag," fragments."
4672 write(iout,*) "Length of each fragment, weight and q0:"
4674 write (iout,*) "Set of restraints #",iset
4676 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4677 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4679 write(iout,*) "constraints between ", npair, "fragments."
4680 write(iout,*) "constraint pairs, weights and q0:"
4682 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4683 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4685 write(iout,*) "angle constraints within ", nfrag_back,&
4686 "backbone fragments."
4687 write(iout,*) "fragment, weights:"
4689 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4690 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4691 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4694 iset=mod(kolor,nset)+1
4697 if(me.eq.king.or..not.out1file) &
4698 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4700 end subroutine read_MDpar
4701 !-----------------------------------------------------------------------------
4705 ! implicit real*8 (a-h,o-z)
4706 ! include 'DIMENSIONS'
4707 ! include 'COMMON.MAP'
4708 ! include 'COMMON.IOUNITS'
4709 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4710 character(len=80) :: mapcard !,ucase
4713 ! real(kind=8) :: var,ene
4716 read (inp,'(a)') mapcard
4717 mapcard=ucase(mapcard)
4718 if (index(mapcard,'PHI').gt.0) then
4720 else if (index(mapcard,'THE').gt.0) then
4722 else if (index(mapcard,'ALP').gt.0) then
4724 else if (index(mapcard,'OME').gt.0) then
4727 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4728 stop 'Error - illegal variable spec in MAP card.'
4730 call readi (mapcard,'RES1',res1(imap),0)
4731 call readi (mapcard,'RES2',res2(imap),0)
4732 if (res1(imap).eq.0) then
4733 res1(imap)=res2(imap)
4734 else if (res2(imap).eq.0) then
4735 res2(imap)=res1(imap)
4737 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4738 write (iout,'(a)') &
4739 'Error - illegal definition of variable group in MAP.'
4740 stop 'Error - illegal definition of variable group in MAP.'
4742 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4743 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4744 call readi(mapcard,'NSTEP',nstep(imap),0)
4745 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4746 write (iout,'(a)') &
4747 'Illegal boundary and/or step size specification in MAP.'
4748 stop 'Illegal boundary and/or step size specification in MAP.'
4752 end subroutine map_read
4753 !-----------------------------------------------------------------------------
4756 use control_data, only: vdisulf
4758 ! implicit real*8 (a-h,o-z)
4759 ! include 'DIMENSIONS'
4760 ! include 'COMMON.IOUNITS'
4761 ! include 'COMMON.GEO'
4762 ! include 'COMMON.CSA'
4763 ! include 'COMMON.BANK'
4764 ! include 'COMMON.CONTROL'
4765 ! character(len=80) :: ucase
4766 character(len=620) :: mcmcard
4768 ! integer :: ntf,ik,iw_pdb
4769 ! real(kind=8) :: var,ene
4771 call card_concat(mcmcard,.true.)
4773 call readi(mcmcard,'NCONF',nconf,50)
4774 call readi(mcmcard,'NADD',nadd,0)
4775 call readi(mcmcard,'JSTART',jstart,1)
4776 call readi(mcmcard,'JEND',jend,1)
4777 call readi(mcmcard,'NSTMAX',nstmax,500000)
4778 call readi(mcmcard,'N0',n0,1)
4779 call readi(mcmcard,'N1',n1,6)
4780 call readi(mcmcard,'N2',n2,4)
4781 call readi(mcmcard,'N3',n3,0)
4782 call readi(mcmcard,'N4',n4,0)
4783 call readi(mcmcard,'N5',n5,0)
4784 call readi(mcmcard,'N6',n6,10)
4785 call readi(mcmcard,'N7',n7,0)
4786 call readi(mcmcard,'N8',n8,0)
4787 call readi(mcmcard,'N9',n9,0)
4788 call readi(mcmcard,'N14',n14,0)
4789 call readi(mcmcard,'N15',n15,0)
4790 call readi(mcmcard,'N16',n16,0)
4791 call readi(mcmcard,'N17',n17,0)
4792 call readi(mcmcard,'N18',n18,0)
4794 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4796 call readi(mcmcard,'NDIFF',ndiff,2)
4797 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4798 call readi(mcmcard,'IS1',is1,1)
4799 call readi(mcmcard,'IS2',is2,8)
4800 call readi(mcmcard,'NRAN0',nran0,4)
4801 call readi(mcmcard,'NRAN1',nran1,2)
4802 call readi(mcmcard,'IRR',irr,1)
4803 call readi(mcmcard,'NSEED',nseed,20)
4804 call readi(mcmcard,'NTOTAL',ntotal,10000)
4805 call reada(mcmcard,'CUT1',cut1,2.0d0)
4806 call reada(mcmcard,'CUT2',cut2,5.0d0)
4807 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4808 call readi(mcmcard,'ICMAX',icmax,3)
4809 call readi(mcmcard,'IRESTART',irestart,0)
4810 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4813 call reada(mcmcard,'DELE',dele,20.0d0)
4814 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4815 call readi(mcmcard,'IREF',iref,0)
4816 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4817 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4818 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4819 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4820 write (iout,*) "NCONF_IN",nconf_in
4822 end subroutine csaread
4823 !-----------------------------------------------------------------------------
4827 use control_data, only: MaxMoveType
4830 ! implicit real*8 (a-h,o-z)
4831 ! include 'DIMENSIONS'
4832 ! include 'COMMON.MCM'
4833 ! include 'COMMON.MCE'
4834 ! include 'COMMON.IOUNITS'
4835 ! character(len=80) :: ucase
4836 character(len=320) :: mcmcard
4839 ! real(kind=8) :: var,ene
4841 call card_concat(mcmcard,.true.)
4842 call readi(mcmcard,'MAXACC',maxacc,100)
4843 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4844 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4845 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4846 call readi(mcmcard,'MAXREPM',maxrepm,200)
4847 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4848 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4849 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4850 call reada(mcmcard,'E_UP',e_up,5.0D0)
4851 call reada(mcmcard,'DELTE',delte,0.1D0)
4852 call readi(mcmcard,'NSWEEP',nsweep,5)
4853 call readi(mcmcard,'NSTEPH',nsteph,0)
4854 call readi(mcmcard,'NSTEPC',nstepc,0)
4855 call reada(mcmcard,'TMIN',tmin,298.0D0)
4856 call reada(mcmcard,'TMAX',tmax,298.0D0)
4857 call readi(mcmcard,'NWINDOW',nwindow,0)
4858 call readi(mcmcard,'PRINT_MC',print_mc,0)
4859 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4860 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4861 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4862 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4863 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4864 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4865 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4866 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4867 if (nwindow.gt.0) then
4868 allocate(winstart(nwindow)) !!el (maxres)
4869 allocate(winend(nwindow)) !!el
4870 allocate(winlen(nwindow)) !!el
4871 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4873 winlen(i)=winend(i)-winstart(i)+1
4876 if (tmax.lt.tmin) tmax=tmin
4877 if (tmax.eq.tmin) then
4881 if (nstepc.gt.0 .and. nsteph.gt.0) then
4882 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4883 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4885 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4886 ! Probabilities of different move types
4887 sumpro_type(0)=0.0D0
4888 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4889 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4890 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4891 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4892 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4893 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4894 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4896 print *,'i',i,' sumprotype',sumpro_type(i)
4897 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4898 print *,'i',i,' sumprotype',sumpro_type(i)
4901 end subroutine mcmread
4902 !-----------------------------------------------------------------------------
4903 subroutine read_minim
4906 ! implicit real*8 (a-h,o-z)
4907 ! include 'DIMENSIONS'
4908 ! include 'COMMON.MINIM'
4909 ! include 'COMMON.IOUNITS'
4910 ! character(len=80) :: ucase
4911 character(len=320) :: minimcard
4913 ! integer :: ntf,ik,iw_pdb
4914 ! real(kind=8) :: var,ene
4916 call card_concat(minimcard,.true.)
4917 call readi(minimcard,'MAXMIN',maxmin,2000)
4918 call readi(minimcard,'MAXFUN',maxfun,5000)
4919 call readi(minimcard,'MINMIN',minmin,maxmin)
4920 call readi(minimcard,'MINFUN',minfun,maxmin)
4921 call reada(minimcard,'TOLF',tolf,1.0D-2)
4922 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4923 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4924 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4925 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4926 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4927 'Options in energy minimization:'
4928 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4929 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4930 'MinMin:',MinMin,' MinFun:',MinFun,&
4931 ' TolF:',TolF,' RTolF:',RTolF
4933 end subroutine read_minim
4934 !-----------------------------------------------------------------------------
4935 subroutine openunits
4937 use MD_data, only: usampl
4940 use control_data, only:out1file
4941 use control, only: getenv_loc
4942 ! implicit real*8 (a-h,o-z)
4943 ! include 'DIMENSIONS'
4946 character(len=16) :: form,nodename
4947 integer :: nodelen,ierror,npos
4949 ! include 'COMMON.SETUP'
4950 ! include 'COMMON.IOUNITS'
4951 ! include 'COMMON.MD'
4952 ! include 'COMMON.CONTROL'
4953 integer :: lenpre,lenpot,lentmp !,ilen
4955 character(len=3) :: out1file_text !,ucase
4956 character(len=3) :: ll
4959 ! integer :: ntf,ik,iw_pdb
4960 ! real(kind=8) :: var,ene
4962 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4963 call getenv_loc("PREFIX",prefix)
4965 call getenv_loc("POT",pot)
4966 call getenv_loc("DIRTMP",tmpdir)
4967 call getenv_loc("CURDIR",curdir)
4968 call getenv_loc("OUT1FILE",out1file_text)
4969 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4970 out1file_text=ucase(out1file_text)
4971 if (out1file_text(1:1).eq."Y") then
4974 out1file=fg_rank.gt.0
4979 if (lentmp.gt.0) then
4980 write (*,'(80(1h!))')
4981 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4982 write (*,'(80(1h!))')
4983 write (*,*)"All output files will be on node /tmp directory."
4985 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4986 if (me.eq.king) then
4987 write (*,*) "The master node is ",nodename
4988 else if (fg_rank.eq.0) then
4989 write (*,*) "I am the CG slave node ",nodename
4991 write (*,*) "I am the FG slave node ",nodename
4994 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4995 lenpre = lentmp+lenpre+1
4997 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4998 ! Get the names and open the input files
4999 #if defined(WINIFL) || defined(WINPGI)
5000 open(1,file=pref_orig(:ilen(pref_orig))// &
5001 '.inp',status='old',readonly,shared)
5002 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5003 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5004 ! Get parameter filenames and open the parameter files.
5005 call getenv_loc('BONDPAR',bondname)
5006 open (ibond,file=bondname,status='old',readonly,shared)
5007 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5008 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5009 call getenv_loc('THETPAR',thetname)
5010 open (ithep,file=thetname,status='old',readonly,shared)
5011 call getenv_loc('ROTPAR',rotname)
5012 open (irotam,file=rotname,status='old',readonly,shared)
5013 call getenv_loc('TORPAR',torname)
5014 open (itorp,file=torname,status='old',readonly,shared)
5015 call getenv_loc('TORDPAR',tordname)
5016 open (itordp,file=tordname,status='old',readonly,shared)
5017 call getenv_loc('FOURIER',fouriername)
5018 open (ifourier,file=fouriername,status='old',readonly,shared)
5019 call getenv_loc('ELEPAR',elename)
5020 open (ielep,file=elename,status='old',readonly,shared)
5021 call getenv_loc('SIDEPAR',sidename)
5022 open (isidep,file=sidename,status='old',readonly,shared)
5024 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5025 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5026 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5027 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5028 call getenv_loc('TORPAR_NUCL',torname_nucl)
5029 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5030 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5031 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5032 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5033 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5036 #elif (defined CRAY) || (defined AIX)
5037 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5039 ! print *,"Processor",myrank," opened file 1"
5040 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5041 ! print *,"Processor",myrank," opened file 9"
5042 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5043 ! Get parameter filenames and open the parameter files.
5044 call getenv_loc('BONDPAR',bondname)
5045 open (ibond,file=bondname,status='old',action='read')
5046 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5047 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5049 ! print *,"Processor",myrank," opened file IBOND"
5050 call getenv_loc('THETPAR',thetname)
5051 open (ithep,file=thetname,status='old',action='read')
5052 ! print *,"Processor",myrank," opened file ITHEP"
5053 call getenv_loc('ROTPAR',rotname)
5054 open (irotam,file=rotname,status='old',action='read')
5055 ! print *,"Processor",myrank," opened file IROTAM"
5056 call getenv_loc('TORPAR',torname)
5057 open (itorp,file=torname,status='old',action='read')
5058 ! print *,"Processor",myrank," opened file ITORP"
5059 call getenv_loc('TORDPAR',tordname)
5060 open (itordp,file=tordname,status='old',action='read')
5061 ! print *,"Processor",myrank," opened file ITORDP"
5062 call getenv_loc('SCCORPAR',sccorname)
5063 open (isccor,file=sccorname,status='old',action='read')
5064 ! print *,"Processor",myrank," opened file ISCCOR"
5065 call getenv_loc('FOURIER',fouriername)
5066 open (ifourier,file=fouriername,status='old',action='read')
5067 ! print *,"Processor",myrank," opened file IFOURIER"
5068 call getenv_loc('ELEPAR',elename)
5069 open (ielep,file=elename,status='old',action='read')
5070 ! print *,"Processor",myrank," opened file IELEP"
5071 call getenv_loc('SIDEPAR',sidename)
5072 open (isidep,file=sidename,status='old',action='read')
5074 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5075 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5076 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5077 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5078 call getenv_loc('TORPAR_NUCL',torname_nucl)
5079 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5080 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5081 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5082 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5083 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5085 call getenv_loc('LIPTRANPAR',liptranname)
5086 open (iliptranpar,file=liptranname,status='old',action='read')
5087 call getenv_loc('TUBEPAR',tubename)
5088 open (itube,file=tubename,status='old',action='read')
5089 call getenv_loc('IONPAR',ionname)
5090 open (iion,file=ionname,status='old',action='read')
5092 ! print *,"Processor",myrank," opened file ISIDEP"
5093 ! print *,"Processor",myrank," opened parameter files"
5095 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5096 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5097 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5098 ! Get parameter filenames and open the parameter files.
5099 call getenv_loc('BONDPAR',bondname)
5100 open (ibond,file=bondname,status='old')
5101 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5102 open (ibond_nucl,file=bondname_nucl,status='old')
5104 call getenv_loc('THETPAR',thetname)
5105 open (ithep,file=thetname,status='old')
5106 call getenv_loc('ROTPAR',rotname)
5107 open (irotam,file=rotname,status='old')
5108 call getenv_loc('TORPAR',torname)
5109 open (itorp,file=torname,status='old')
5110 call getenv_loc('TORDPAR',tordname)
5111 open (itordp,file=tordname,status='old')
5112 call getenv_loc('SCCORPAR',sccorname)
5113 open (isccor,file=sccorname,status='old')
5114 call getenv_loc('FOURIER',fouriername)
5115 open (ifourier,file=fouriername,status='old')
5116 call getenv_loc('ELEPAR',elename)
5117 open (ielep,file=elename,status='old')
5118 call getenv_loc('SIDEPAR',sidename)
5119 open (isidep,file=sidename,status='old')
5121 open (ithep_nucl,file=thetname_nucl,status='old')
5122 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5123 open (irotam_nucl,file=rotname_nucl,status='old')
5124 call getenv_loc('TORPAR_NUCL',torname_nucl)
5125 open (itorp_nucl,file=torname_nucl,status='old')
5126 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5127 open (itordp_nucl,file=tordname_nucl,status='old')
5128 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5129 open (isidep_nucl,file=sidename_nucl,status='old')
5131 call getenv_loc('LIPTRANPAR',liptranname)
5132 open (iliptranpar,file=liptranname,status='old')
5133 call getenv_loc('TUBEPAR',tubename)
5134 open (itube,file=tubename,status='old')
5135 call getenv_loc('IONPAR',ionname)
5136 open (iion,file=ionname,status='old')
5138 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5140 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5141 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5142 ! Get parameter filenames and open the parameter files.
5143 call getenv_loc('BONDPAR',bondname)
5144 open (ibond,file=bondname,status='old',action='read')
5145 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5146 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5147 call getenv_loc('THETPAR',thetname)
5148 open (ithep,file=thetname,status='old',action='read')
5149 call getenv_loc('ROTPAR',rotname)
5150 open (irotam,file=rotname,status='old',action='read')
5151 call getenv_loc('TORPAR',torname)
5152 open (itorp,file=torname,status='old',action='read')
5153 call getenv_loc('TORDPAR',tordname)
5154 open (itordp,file=tordname,status='old',action='read')
5155 call getenv_loc('SCCORPAR',sccorname)
5156 open (isccor,file=sccorname,status='old',action='read')
5158 call getenv_loc('THETPARPDB',thetname_pdb)
5159 print *,"thetname_pdb ",thetname_pdb
5160 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5161 print *,ithep_pdb," opened"
5163 call getenv_loc('FOURIER',fouriername)
5164 open (ifourier,file=fouriername,status='old',readonly)
5165 call getenv_loc('ELEPAR',elename)
5166 open (ielep,file=elename,status='old',readonly)
5167 call getenv_loc('SIDEPAR',sidename)
5168 open (isidep,file=sidename,status='old',readonly)
5170 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5171 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5172 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5173 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5174 call getenv_loc('TORPAR_NUCL',torname_nucl)
5175 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5176 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5177 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5178 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5179 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5180 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5181 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5182 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5183 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5184 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5185 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5186 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5187 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5190 call getenv_loc('LIPTRANPAR',liptranname)
5191 open (iliptranpar,file=liptranname,status='old',action='read')
5192 call getenv_loc('TUBEPAR',tubename)
5193 open (itube,file=tubename,status='old',action='read')
5194 call getenv_loc('IONPAR',ionname)
5195 open (iion,file=ionname,status='old',action='read')
5198 call getenv_loc('ROTPARPDB',rotname_pdb)
5199 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5202 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5203 #if defined(WINIFL) || defined(WINPGI)
5204 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5205 #elif (defined CRAY) || (defined AIX)
5206 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5208 open (iscpp_nucl,file=scpname_nucl,status='old')
5210 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5215 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5216 ! Use -DOLDSCP to use hard-coded constants instead.
5218 call getenv_loc('SCPPAR',scpname)
5219 #if defined(WINIFL) || defined(WINPGI)
5220 open (iscpp,file=scpname,status='old',readonly,shared)
5221 #elif (defined CRAY) || (defined AIX)
5222 open (iscpp,file=scpname,status='old',action='read')
5224 open (iscpp,file=scpname,status='old')
5226 open (iscpp,file=scpname,status='old',action='read')
5229 call getenv_loc('PATTERN',patname)
5230 #if defined(WINIFL) || defined(WINPGI)
5231 open (icbase,file=patname,status='old',readonly,shared)
5232 #elif (defined CRAY) || (defined AIX)
5233 open (icbase,file=patname,status='old',action='read')
5235 open (icbase,file=patname,status='old')
5237 open (icbase,file=patname,status='old',action='read')
5240 ! Open output file only for CG processes
5241 ! print *,"Processor",myrank," fg_rank",fg_rank
5242 if (fg_rank.eq.0) then
5244 if (nodes.eq.1) then
5247 npos = dlog10(dfloat(nodes-1))+1
5249 if (npos.lt.3) npos=3
5250 write (liczba,'(i1)') npos
5251 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5253 write (liczba,form) me
5254 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5255 liczba(:ilen(liczba))
5256 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5258 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5260 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5261 liczba(:ilen(liczba))//'.mol2'
5262 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5263 liczba(:ilen(liczba))//'.stat'
5265 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5266 //liczba(:ilen(liczba))//'.stat')
5267 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5270 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5271 liczba(:ilen(liczba))//'.const'
5276 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5277 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5278 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5279 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5280 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5282 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5284 rest2name=prefix(:ilen(prefix))//'.rst'
5286 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5289 #if defined(AIX) || defined(PGI)
5290 if (me.eq.king .or. .not. out1file) &
5291 open(iout,file=outname,status='unknown')
5293 if (fg_rank.gt.0) then
5294 write (liczba,'(i3.3)') myrank/nfgtasks
5295 write (ll,'(bz,i3.3)') fg_rank
5296 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5301 open(igeom,file=intname,status='unknown',position='append')
5302 open(ipdb,file=pdbname,status='unknown')
5303 open(imol2,file=mol2name,status='unknown')
5304 open(istat,file=statname,status='unknown',position='append')
5306 !1out open(iout,file=outname,status='unknown')
5309 if (me.eq.king .or. .not.out1file) &
5310 open(iout,file=outname,status='unknown')
5312 if (fg_rank.gt.0) then
5313 write (liczba,'(i3.3)') myrank/nfgtasks
5314 write (ll,'(bz,i3.3)') fg_rank
5315 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5320 open(igeom,file=intname,status='unknown',access='append')
5321 open(ipdb,file=pdbname,status='unknown')
5322 open(imol2,file=mol2name,status='unknown')
5323 open(istat,file=statname,status='unknown',access='append')
5325 !1out open(iout,file=outname,status='unknown')
5328 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5329 csa_seed=prefix(:lenpre)//'.CSA.seed'
5330 csa_history=prefix(:lenpre)//'.CSA.history'
5331 csa_bank=prefix(:lenpre)//'.CSA.bank'
5332 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5333 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5334 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5335 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5336 csa_int=prefix(:lenpre)//'.int'
5337 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5338 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5339 csa_in=prefix(:lenpre)//'.CSA.in'
5340 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5343 write (iout,'(80(1h-))')
5344 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5345 write (iout,'(80(1h-))')
5346 write (iout,*) "Input file : ",&
5347 pref_orig(:ilen(pref_orig))//'.inp'
5348 write (iout,*) "Output file : ",&
5349 outname(:ilen(outname))
5351 write (iout,*) "Sidechain potential file : ",&
5352 sidename(:ilen(sidename))
5354 write (iout,*) "SCp potential file : ",&
5355 scpname(:ilen(scpname))
5357 write (iout,*) "Electrostatic potential file : ",&
5358 elename(:ilen(elename))
5359 write (iout,*) "Cumulant coefficient file : ",&
5360 fouriername(:ilen(fouriername))
5361 write (iout,*) "Torsional parameter file : ",&
5362 torname(:ilen(torname))
5363 write (iout,*) "Double torsional parameter file : ",&
5364 tordname(:ilen(tordname))
5365 write (iout,*) "SCCOR parameter file : ",&
5366 sccorname(:ilen(sccorname))
5367 write (iout,*) "Bond & inertia constant file : ",&
5368 bondname(:ilen(bondname))
5369 write (iout,*) "Bending parameter file : ",&
5370 thetname(:ilen(thetname))
5371 write (iout,*) "Rotamer parameter file : ",&
5372 rotname(:ilen(rotname))
5375 write (iout,*) "Thetpdb parameter file : ",&
5376 thetname_pdb(:ilen(thetname_pdb))
5379 write (iout,*) "Threading database : ",&
5380 patname(:ilen(patname))
5382 write (iout,*)" DIRTMP : ",&
5384 write (iout,'(80(1h-))')
5387 end subroutine openunits
5388 !-----------------------------------------------------------------------------
5391 use geometry_data, only: nres,dc
5393 ! implicit real*8 (a-h,o-z)
5394 ! include 'DIMENSIONS'
5395 ! include 'COMMON.CHAIN'
5396 ! include 'COMMON.IOUNITS'
5397 ! include 'COMMON.MD'
5400 ! real(kind=8) :: var,ene
5402 open(irest2,file=rest2name,status='unknown')
5403 read(irest2,*) totT,EK,potE,totE,t_bath
5406 ! AL 4/17/17: Now reading d_t(0,:) too
5408 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5411 ! AL 4/17/17: Now reading d_c(0,:) too
5413 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5416 read (irest2,*) iset
5420 end subroutine readrst
5421 !-----------------------------------------------------------------------------
5422 subroutine read_fragments
5426 use control_data, only:out1file
5429 ! implicit real*8 (a-h,o-z)
5430 ! include 'DIMENSIONS'
5434 ! include 'COMMON.SETUP'
5435 ! include 'COMMON.CHAIN'
5436 ! include 'COMMON.IOUNITS'
5437 ! include 'COMMON.MD'
5438 ! include 'COMMON.CONTROL'
5441 ! real(kind=8) :: var,ene
5443 read(inp,*) nset,nfrag,npair,nfrag_back
5445 !el from module energy
5446 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5447 if(.not.allocated(wfrag_back)) then
5448 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5449 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5451 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5452 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5454 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5455 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5458 if(me.eq.king.or..not.out1file) &
5459 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5460 " nfrag_back",nfrag_back
5462 read(inp,*) mset(iset)
5464 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5466 if(me.eq.king.or..not.out1file) &
5467 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5468 ifrag(2,i,iset), qinfrag(i,iset)
5471 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5473 if(me.eq.king.or..not.out1file) &
5474 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5475 ipair(2,i,iset), qinpair(i,iset)
5478 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5479 wfrag_back(3,i,iset),&
5480 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5481 if(me.eq.king.or..not.out1file) &
5482 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5483 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5487 end subroutine read_fragments
5488 !-----------------------------------------------------------------------------
5490 !-----------------------------------------------------------------------------
5494 ! implicit real*8 (a-h,o-z)
5495 ! include 'DIMENSIONS'
5496 ! include 'COMMON.CSA'
5497 ! include 'COMMON.BANK'
5498 ! include 'COMMON.IOUNITS'
5500 ! integer :: ntf,ik,iw_pdb
5501 ! real(kind=8) :: var,ene
5503 open(icsa_in,file=csa_in,status="old",err=100)
5504 read(icsa_in,*) nconf
5505 read(icsa_in,*) jstart,jend
5506 read(icsa_in,*) nstmax
5507 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5508 read(icsa_in,*) nran0,nran1,irr
5509 read(icsa_in,*) nseed
5510 read(icsa_in,*) ntotal,cut1,cut2
5511 read(icsa_in,*) estop
5512 read(icsa_in,*) icmax,irestart
5513 read(icsa_in,*) ntbankm,dele,difcut
5514 read(icsa_in,*) iref,rmscut,pnccut
5515 read(icsa_in,*) ndiff
5522 end subroutine csa_read
5523 !-----------------------------------------------------------------------------
5524 subroutine initial_write
5527 ! implicit real*8 (a-h,o-z)
5528 ! include 'DIMENSIONS'
5529 ! include 'COMMON.CSA'
5530 ! include 'COMMON.BANK'
5531 ! include 'COMMON.IOUNITS'
5533 ! integer :: ntf,ik,iw_pdb
5534 ! real(kind=8) :: var,ene
5536 open(icsa_seed,file=csa_seed,status="unknown")
5537 write(icsa_seed,*) "seed"
5539 #if defined(AIX) || defined(PGI)
5540 open(icsa_history,file=csa_history,status="unknown",&
5543 open(icsa_history,file=csa_history,status="unknown",&
5546 write(icsa_history,*) nconf
5547 write(icsa_history,*) jstart,jend
5548 write(icsa_history,*) nstmax
5549 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5550 write(icsa_history,*) nran0,nran1,irr
5551 write(icsa_history,*) nseed
5552 write(icsa_history,*) ntotal,cut1,cut2
5553 write(icsa_history,*) estop
5554 write(icsa_history,*) icmax,irestart
5555 write(icsa_history,*) ntbankm,dele,difcut
5556 write(icsa_history,*) iref,rmscut,pnccut
5557 write(icsa_history,*) ndiff
5559 write(icsa_history,*)
5562 open(icsa_bank1,file=csa_bank1,status="unknown")
5563 write(icsa_bank1,*) 0
5567 end subroutine initial_write
5568 !-----------------------------------------------------------------------------
5569 subroutine restart_write
5572 ! implicit real*8 (a-h,o-z)
5573 ! include 'DIMENSIONS'
5574 ! include 'COMMON.IOUNITS'
5575 ! include 'COMMON.CSA'
5576 ! include 'COMMON.BANK'
5578 ! integer :: ntf,ik,iw_pdb
5579 ! real(kind=8) :: var,ene
5581 #if defined(AIX) || defined(PGI)
5582 open(icsa_history,file=csa_history,position="append")
5584 open(icsa_history,file=csa_history,access="append")
5586 write(icsa_history,*)
5587 write(icsa_history,*) "This is restart"
5588 write(icsa_history,*)
5589 write(icsa_history,*) nconf
5590 write(icsa_history,*) jstart,jend
5591 write(icsa_history,*) nstmax
5592 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5593 write(icsa_history,*) nran0,nran1,irr
5594 write(icsa_history,*) nseed
5595 write(icsa_history,*) ntotal,cut1,cut2
5596 write(icsa_history,*) estop
5597 write(icsa_history,*) icmax,irestart
5598 write(icsa_history,*) ntbankm,dele,difcut
5599 write(icsa_history,*) iref,rmscut,pnccut
5600 write(icsa_history,*) ndiff
5601 write(icsa_history,*)
5602 write(icsa_history,*) "irestart is: ", irestart
5604 write(icsa_history,*)
5608 end subroutine restart_write
5609 !-----------------------------------------------------------------------------
5611 !-----------------------------------------------------------------------------
5612 subroutine write_pdb(npdb,titelloc,ee)
5614 ! implicit real*8 (a-h,o-z)
5615 ! include 'DIMENSIONS'
5616 ! include 'COMMON.IOUNITS'
5617 character(len=50) :: titelloc1
5618 character*(*) :: titelloc
5619 character(len=3) :: zahl
5620 character(len=5) :: liczba5
5622 integer :: npdb !,ilen
5626 ! real(kind=8) :: var,ene
5630 if (npdb.lt.1000) then
5631 call numstr(npdb,zahl)
5632 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5634 if (npdb.lt.10000) then
5635 write(liczba5,'(i1,i4)') 0,npdb
5637 write(liczba5,'(i5)') npdb
5639 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5641 call pdbout(ee,titelloc1,ipdb)
5644 end subroutine write_pdb
5645 !-----------------------------------------------------------------------------
5647 !-----------------------------------------------------------------------------
5648 subroutine write_thread_summary
5649 ! Thread the sequence through a database of known structures
5650 use control_data, only: refstr
5652 use energy_data, only: n_ene_comp
5654 ! implicit real*8 (a-h,o-z)
5655 ! include 'DIMENSIONS'
5657 use MPI_data !include 'COMMON.INFO'
5660 ! include 'COMMON.CONTROL'
5661 ! include 'COMMON.CHAIN'
5662 ! include 'COMMON.DBASE'
5663 ! include 'COMMON.INTERACT'
5664 ! include 'COMMON.VAR'
5665 ! include 'COMMON.THREAD'
5666 ! include 'COMMON.FFIELD'
5667 ! include 'COMMON.SBRIDGE'
5668 ! include 'COMMON.HEADER'
5669 ! include 'COMMON.NAMES'
5670 ! include 'COMMON.IOUNITS'
5671 ! include 'COMMON.TIME1'
5673 integer,dimension(maxthread) :: ip
5674 real(kind=8),dimension(0:n_ene) :: energia
5676 integer :: i,j,ii,jj,ipj,ik,kk,ist
5677 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5679 write (iout,'(30x,a/)') &
5680 ' *********** Summary threading statistics ************'
5681 write (iout,'(a)') 'Initial energies:'
5682 write (iout,'(a4,2x,a12,14a14,3a8)') &
5683 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5684 'RMSnat','NatCONT','NNCONT','RMS'
5685 ! Energy sort patterns
5690 enet=ener(n_ene-1,ip(i))
5693 if (ener(n_ene-1,ip(j)).lt.enet) then
5695 enet=ener(n_ene-1,ip(j))
5707 ist=nres_base(2,ii)+ipatt(2,i)
5709 energia(i)=ener0(kk,i)
5711 etot=ener0(n_ene_comp+1,i)
5712 rmsnat=ener0(n_ene_comp+2,i)
5713 rms=ener0(n_ene_comp+3,i)
5714 frac=ener0(n_ene_comp+4,i)
5715 frac_nn=ener0(n_ene_comp+5,i)
5718 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5719 i,str_nam(ii),ist+1,&
5720 (energia(print_order(kk)),kk=1,nprint_ene),&
5721 etot,rmsnat,frac,frac_nn,rms
5723 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5724 i,str_nam(ii),ist+1,&
5725 (energia(print_order(kk)),kk=1,nprint_ene),etot
5728 write (iout,'(//a)') 'Final energies:'
5729 write (iout,'(a4,2x,a12,17a14,3a8)') &
5730 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5731 'RMSnat','NatCONT','NNCONT','RMS'
5735 ist=nres_base(2,ii)+ipatt(2,i)
5737 energia(kk)=ener(kk,ik)
5739 etot=ener(n_ene_comp+1,i)
5740 rmsnat=ener(n_ene_comp+2,i)
5741 rms=ener(n_ene_comp+3,i)
5742 frac=ener(n_ene_comp+4,i)
5743 frac_nn=ener(n_ene_comp+5,i)
5744 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5745 i,str_nam(ii),ist+1,&
5746 (energia(print_order(kk)),kk=1,nprint_ene),&
5747 etot,rmsnat,frac,frac_nn,rms
5749 write (iout,'(/a/)') 'IEXAM array:'
5750 write (iout,'(i5)') nexcl
5752 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5754 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5755 'Max. time for threading step ',max_time_for_thread,&
5756 'Average time for threading step: ',ave_time_for_thread
5758 end subroutine write_thread_summary
5759 !-----------------------------------------------------------------------------
5760 subroutine write_stat_thread(ithread,ipattern,ist)
5762 use energy_data, only: n_ene_comp
5764 ! implicit real*8 (a-h,o-z)
5765 ! include "DIMENSIONS"
5766 ! include "COMMON.CONTROL"
5767 ! include "COMMON.IOUNITS"
5768 ! include "COMMON.THREAD"
5769 ! include "COMMON.FFIELD"
5770 ! include "COMMON.DBASE"
5771 ! include "COMMON.NAMES"
5772 real(kind=8),dimension(0:n_ene) :: energia
5774 integer :: ithread,ipattern,ist,i
5775 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5777 #if defined(AIX) || defined(PGI)
5778 open(istat,file=statname,position='append')
5780 open(istat,file=statname,access='append')
5783 energia(i)=ener(i,ithread)
5785 etot=ener(n_ene_comp+1,ithread)
5786 rmsnat=ener(n_ene_comp+2,ithread)
5787 rms=ener(n_ene_comp+3,ithread)
5788 frac=ener(n_ene_comp+4,ithread)
5789 frac_nn=ener(n_ene_comp+5,ithread)
5790 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5791 ithread,str_nam(ipattern),ist+1,&
5792 (energia(print_order(i)),i=1,nprint_ene),&
5793 etot,rmsnat,frac,frac_nn,rms
5796 end subroutine write_stat_thread
5797 !-----------------------------------------------------------------------------
5799 !-----------------------------------------------------------------------------
5800 end module io_config